Blaise Delmotte
email: blaise.delmotte@imft.fr

Eric Climent
email: eric.climent@imft.fr

Franck Plouraboué
email: franck.plouraboue@imft.fr

A general formulation of Bead Models applied to flexible fibers and active filaments at low Reynolds number

Keywords:

This contribution provides a general framework to use Lagrange multipliers for the simulation of low Reynolds number fiber dynamics based on Bead Models (BM). This formalism provides an efficient method to account for kinematic constraints. We illustrate, with several examples, to which extent the proposed formulation offers a flexible and versatile framework for the quantitative modeling of flexible fibers deformation and rotation in shear flow, the dynamics of actuated filaments and the propulsion of active swimmers. Furthermore, a new contact model called Gears Model is proposed and successfully tested. It avoids the use of numerical artifices such as repulsive forces between adjacent beads, a source of numerical difficulties in the temporal integration of previous Bead Models.

Introduction

The dynamics of solid-liquid suspensions is a longstanding topic of research while it combines difficulties arising from the coupling of multi-body interactions in a viscous fluid with possible deformations of flexible objects such as fibers. A vast literature exists on the response of suspensions of solid spherical or non-spherical particles due to its ubiquitous interest in natural and industrial processes. When the objects have the ability to deform many complications arise. The coupling between suspended particles will depend on the positions (possibly orientations) but also on the shape of individuals, introducing intricate effects of the history of the suspension.

When the aspect ratio of deformable objects is large, those are generally designated as fibers. Many previous investigations of fiber dynamics, have focused on the dynamics of rigid fibers or rods [START_REF] Cox | The motion of long slender bodies in a viscous fluid. Part 2. Shear flow[END_REF][START_REF] Meunier | Friction coefficient of rod-like chains of spheres at very low Reynolds numbers. II. Numerical simulations[END_REF]. Compared to the very large number of references related to particle suspensions, lower attention has been paid to the more complicated system of flexible fibers in a fluid.

Suspension of flexible fibers are encountered in the study of polymer dynamics [START_REF] Yamakawa | Transport properties of polymer chains in dilute solution: hydrodynamic interaction[END_REF][START_REF] Yamanoi | Analysis of rheological properties of fibre suspensions in a Newtonian fluid by direct fibre simulation. Part 2: Flexible fibre suspensions[END_REF] whose rheology depends on the formation of networks and the occurrence of entanglement. The motion of fibers in a viscous fluid has a strong effect on its bulk viscosity, microstructure, drainage rate, filtration ability, and flocculation properties. The dynamic response of such complex solutions is still an open issue while time-dependent structural changes of the dispersed fibers can dramatically modify the overall process (such as operation units in wood pulp and paper industry, flow molding techniques of composites, water purification). Biological fibers such as DNA or actin filaments have also attracted many researches to understand the relation between flexibility and physiological properties [START_REF] Jian | A combined Wormlike-Chain and Bead Model for dynamic simulations of long linear DNA[END_REF].

Flexible fibers do not only passively respond to carrying flow gradients but can also be dynamically activated. Many of single cell micro-organisms that propel themselves in a fluid utilize a long flagellum tail connected to the cell body. Spermatozoa (and more generally one-armed swimmers) swim by propagating bending waves along their flagellum tail to generate a net translation using cyclic non-reciprocal strategy at low Reynolds number [START_REF] Purcell | Life at low Reynolds number[END_REF]. These natural swimmers have been modeled by artificial swimmers (joint microbeads) actuated by an oscillating ambient electric or magnetic field which opens breakthrough technologies for drug on-demand delivery in the human body [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF].

Many numerical methods have been proposed to tackle elasto-hydrodynamic coupling between a fluid flow and deformable objects, i.e. the balance between viscous drag and elastic stresses. Among those, "mesh-oriented" approaches have the ambition of solving a complete continuum mechanics description of the fluid/solid interaction, even though some approximations are mandatory to describe those at the fluid/solid interface. Without being all-comprehensive, one can cite immerse boundary methods (e.g. [START_REF] Fogelson | A fast numerical method for solving the three-dimensional Stokes' equations in the presence of suspended particles[END_REF][START_REF] Stockie | Simulating the motion of flexible pulp fibres using the immersed boundary method[END_REF][START_REF] Zhu | Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[END_REF][START_REF] Zhu | Drag of a flexible fiber in a 2D moving viscous fluid[END_REF]), extended finite elements (e.g. [START_REF] Wagner | The extended finite element method for rigid particles in Stokes flow[END_REF]), penalty methods [START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF][START_REF] Decoene | Microscopic modelling of active bacterial suspensions[END_REF], particle-mesh Ewald methods [START_REF] Saintillan | A smooth particle-mesh Ewald algorithm for Stokes suspension simulations: the sedimentation of fibers[END_REF], regularized Stokeslets [START_REF] Olson | Modeling the dynamics of an elastic rod with intrinsic curvature and twist using a regularized Stokes formulation[END_REF][START_REF] Simons | The dynamics of sperm detachment from epithelium in a coupled fluid-biochemical model of hyperactivated motility[END_REF], Force Coupling Method [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF].

In the specific context of low Reynolds number elastohydrodynamics [START_REF] Wiggins | Flexive and propulsive dynamics of elastica at low Reynolds number[END_REF], difficulties arise when numerically solving the dynamics of rigid objects since the time scale associated with elastic waves propagation within the solid can be similar to the viscous dissipation time-scale. In the context of self propelled objects the ratio of these time scales is called "Sperm number". When the Sperm number is smaller or equal to one, the object temporal response is stiff, and requires small time steps to capture fast deformation modes. In this regime, fluid/structure interaction effects are difficult to capture. One possible way to circumvent such difficulties is to use the knowledge of hydrodynamic interactions of simple objects in Stokes flow.

This strategy is the one pursued by the Bead Model (BM) whose aim is to describe a complex deformable object by the flexible assembly of simple rigid ones. Such flexible assemblies are generally composed of beads (spheres or ellipsoids) interacting by some elastic and repulsive forces, as well as with the surrounding fluid. For long elongated structures, alternative approaches to BM are indeed possible such as slender body approximation [START_REF] Cox | The motion of long slender bodies in a viscous fluid. Part 2. Shear flow[END_REF][START_REF] Tornberg | Simulating the dynamics and interactions of flexible fibers in Stokes flows[END_REF][START_REF] Gaffney | Mammalian sperm motility: observation and theory[END_REF][START_REF] Pozrikidis | Shear flow past slender elastic rods attached to a plane[END_REF] or Resistive Force Theory [START_REF] Lauga | Floppy swimming: viscous locomotion of actuated elastica[END_REF][START_REF] Coq | Helical beating of an actuated elastic filament[END_REF][START_REF] Gadêlha | Nonlinear instability in flagellar dynamics: a novel modulation mechanism in sperm migration?[END_REF].

One important advantage of BM which might explain their use among various communities (polymer Physics [START_REF] Meunier | Friction coefficient of rod-like chains of spheres at very low Reynolds numbers. II. Numerical simulations[END_REF][START_REF] Yamakawa | Transport properties of polymer chains in dilute solution: hydrodynamic interaction[END_REF][START_REF] Yamanoi | Analysis of rheological properties of fibre suspensions in a Newtonian fluid by direct fibre simulation. Part 2: Flexible fibre suspensions[END_REF][START_REF] Jian | A combined Wormlike-Chain and Bead Model for dynamic simulations of long linear DNA[END_REF][START_REF] Gao | Numerical and experimental study of a rotating magnetic particle chain in a viscous fluid[END_REF][START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF][START_REF] Jendrejack | Hydrodynamic interactions in long chain polymers: application of the Chebyshev polynomial approximation in stochastic simulations[END_REF][START_REF] Jendrejack | Stochastic simulations of DNA in flow: dynamics and the effects of hydrodynamic interactions[END_REF][START_REF] Schroeder | Effect of hydrodynamic interactions on DNA dynamics in extensional flow: simulation and single molecule experiment[END_REF][START_REF] Montesi | Brownian dynamics algorithm for bead-rod semiflexible chain with anisotropic friction[END_REF][START_REF] Schlagberger | Orientation of elastic rods in homogeneous Stokes flow[END_REF][START_REF] Bailey | Efficient constraint dynamics using MILC SHAKE[END_REF][START_REF] Yamanoi | Stokesian dynamics simulation of the role of hydrodynamic interactions on the behavior of a single particle suspending in a Newtonian fluid. Part 1. 1D flexible and rigid fibers[END_REF], micro-swimmer modeling in bio-fluid mechanics [START_REF] Wada | Non-equilibrium hydrodynamics of a rotating filament[END_REF][START_REF] Manghi | Propulsion with a rotating elastic nanorod[END_REF][START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF][START_REF] Swan | Modeling hydrodynamic self-propulsion with Stokesian Dynamics. Or teaching Stokesian Dynamics to swim[END_REF][START_REF] Bilbao | Navigation and chemotaxis of nematodes in bulk and confined fluids[END_REF][START_REF] Gauger | Numerical study of a microscopic artificial swimmer[END_REF][START_REF] Keaveny | Dynamics of structures in active suspensions of paramagnetic particles and applications to artificial micro-swimmers[END_REF][START_REF] Lowe | Dynamics of filaments: modelling the dynamics of driven microfilaments[END_REF][START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF][START_REF] Berman | Undulatory locomotion of finite filaments: lessons from Caenorhabditis elegans[END_REF], flexible fiber in chemical engineering [START_REF] Joung | Direct simulation of flexible fibers[END_REF][START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF][START_REF] Qi | Direct simulations of flexible cylindrical fiber suspensions in finite Reynolds number flows[END_REF][START_REF] Ross | Dynamic simulation of flexible fibers composed of linked rigid bodies[END_REF][START_REF] Skjetne | Simulation of single fiber dynamics[END_REF][START_REF] Slowicka | Dynamics of fibers in a wide microchannel[END_REF][START_REF] Hoffman | The dynamics of the coil-stretch transition for long, flexible polymers in planar mixed flows[END_REF][START_REF] Wang | Flipping, scooping, and spinning: drift of rigid curved nonchiral fibers in simple shear flow[END_REF]), relies on their parametric versatility, their ubiquitous character and their relative easy implementation. We provide a deeper, comparative and critical discussion about BM in Section 2. However, we would like to stress that the presented model is more clearly oriented toward micro-swimmer modeling than polymer dynamics.

One should also add that BM can be coupled to mesh-oriented approaches in order to provide accurate description of hydrodynamic interactions among large collection of deformable objects at moderate numerical cost [START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF]. Many authors only consider free drain, i.e. no Hydrodynamic Interactions (no HI), [START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF][START_REF] Ross | Dynamic simulation of flexible fibers composed of linked rigid bodies[END_REF][START_REF] Skjetne | Simulation of single fiber dynamics[END_REF][START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF] or far field interactions associated with the Rotne-Prager-Yamakawa tensor [START_REF] Wada | Non-equilibrium hydrodynamics of a rotating filament[END_REF][START_REF] Manghi | Propulsion with a rotating elastic nanorod[END_REF][START_REF] Gauger | Numerical study of a microscopic artificial swimmer[END_REF][START_REF] Wajnryb | Generalization of the Rotne-Prager-Yamakawa mobility and shear disturbance tensors[END_REF]. This is supported by the fact that far-field hydrodynamic interactions already provide accurate predictions for the dynamics of a single flexible fiber when compared to experimental observations or numerical results. In order to illustrate the method we use, for convenience, the Rotne-Prager-Yamakawa tensor to model hydrodynamic interactions. We wish to stress here that this is not a limitation of the presented method, since the presented formulation holds for any mobility problem formulation. However, it turns out that for each configuration we tested, our model gave very good comparisons with other predictions, including those providing more accurate description of the hydrodynamic interactions.

The paper is organized as follows. First, we give a detailed presentation of the Bead Model for the simulation of flexible fibers. In this section, we propose a general formulation of kinematic constraints using the framework of Lagrange multipliers. This general formulation is used to present a new Bead Model, namely the Gears Model which surpasses existing models on numerical aspects. The second part of the paper is devoted to comparisons and validations of Bead Models for different configurations of flexible fibers (experiencing a flow or actuated filaments).

Finally, we conclude the paper by summarizing the achievements we obtain with our model and open new perspectives to this work.

The Bead Model

Detailed review of previous Bead Models

The Bead Model (BM) aims at discretizing any flexible object with interacting beads. Interactions between beads break down into three categories: hydrodynamic interactions, elastic and kinematic constraint forces. Hydrodynamics of the whole object result from multibody hydrodynamic interactions between beads. In the context of low Reynolds number, the relationship between stresses and velocities is linear. Thus, the velocity of the assembly depends linearly on the forces and torques applied on each of its elements. Elastic forces and torques are prescribed according to classical elasticity theory [START_REF] Landau | Theory of Elasticity[END_REF] of flexible matter. Constraint forces ensure that the beads obey any imposed kinematic constraint, e.g. fixed distance between adjacent particles. All of these interactions can be treated separately as long as they are addressed in a consistent order. The latter is the cornerstone which differentiates previous works in the literature from ours. Numerous strategies have been employed to handle kinematic constraints. [START_REF] Schlagberger | Orientation of elastic rods in homogeneous Stokes flow[END_REF][START_REF] Yamanoi | Stokesian dynamics simulation of the role of hydrodynamic interactions on the behavior of a single particle suspending in a Newtonian fluid. Part 1. 1D flexible and rigid fibers[END_REF][START_REF] Wada | Non-equilibrium hydrodynamics of a rotating filament[END_REF][START_REF] Gauger | Numerical study of a microscopic artificial swimmer[END_REF] and [START_REF] Slowicka | Dynamics of fibers in a wide microchannel[END_REF] used a linear spring to model the resistance to stretching and compression without any constraint on the bead rotational motion (Fig. 1). The resulting stretching force reads:

F s = -k s r i,i+1 -r 0 i,i+1 (1)
where • k s is the spring stiffness, • r i,i+1 = r i+1r i is the distance vector between two adjacent beads (for simplicity, equations and figures will be presented for beads 1 and 2 and can easily be generalized to beads i and i + 1), • r 0 1,2 is the vector corresponding to equilibrium. However, regarding the connectivity constraint, the spring model is somehow approximate. A linear spring is prone to uncontrolled oscillations and the problem may become unstable. Many other authors, among which [START_REF] Jendrejack | Hydrodynamic interactions in long chain polymers: application of the Chebyshev polynomial approximation in stochastic simulations[END_REF][START_REF] Jendrejack | Stochastic simulations of DNA in flow: dynamics and the effects of hydrodynamic interactions[END_REF][START_REF] Schroeder | Effect of hydrodynamic interactions on DNA dynamics in extensional flow: simulation and single molecule experiment[END_REF], thus use non-linear spring models for a better description of polymer physics. Nevertheless, the repulsive force stiffness has an important numerical cost in time-stepping as will be discussed in Section 2.6.3. Furthermore, unconstrained bead rotational motion leads to spurious hydrodynamic interactions and thus limits the range of applications for these BM.

Alternatively, [START_REF] Qi | Direct simulations of flexible cylindrical fiber suspensions in finite Reynolds number flows[END_REF][START_REF] Ross | Dynamic simulation of flexible fibers composed of linked rigid bodies[END_REF][START_REF] Skjetne | Simulation of single fiber dynamics[END_REF][START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF] and [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF] constrained the motion of the beads such that the contact point for each pair c i remains the same. While more representative of a flexible object, this approach exhibits two main drawbacks:

1. a gap between beads is necessary to allow the object to bend (see Fig. 2), 2. it requires an additional center to center repulsive force, and thus more tuning numerical parameters to prevent overlapping between adjacent beads.

Consider two adjacent beads, with radius a, linked by a hinge c 1 (typically called ball and socket joint). The gap ε g defines the distance between the sphere surfaces and the joint (see Fig. 3). Denote p i the vector attached to bead i pointing towards the next joint, i.e. the contact point c i . The connectivity between two contiguous bodies writes:

r 1 + (a + ε g )p 1 -r 2 -(a + ε g )p 2 = 0 (2)
and its time derivative

ṙ1 -(a + ε g )p 1 × ω 1 -ṙ2 + (a + ε g )p 2 × ω 2 = 0. (3)
ṙi and ω i are the translational and rotational velocities of bead i. The constraint forces and torques associated to (3) are obtained either by solving a linear system of equations involving beads velocities [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF], or by inserting (3) into the equations of motion when neglecting hydrodynamic interactions [START_REF] Ross | Dynamic simulation of flexible fibers composed of linked rigid bodies[END_REF][START_REF] Skjetne | Simulation of single fiber dynamics[END_REF].

The gap width 2ε g controls the maximum curvature κ J max allowed without overlapping. From the sine rule, one can derive the simple equation relating ε g and κ J max

κ J max = 1 -( a a+ε g ) 2 a ( 4
)
Once aware of these limitations, the gap ε g , range and strength of the repulsive force should be prescribed depending on the problem to be addressed. [START_REF] Keaveny | Spiral swimming of an artificial micro-swimmer[END_REF] and [START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF] proposed a more sophisticated Joint Model than those hitherto cited, using a full description of the links dynamics along the curvilinear abscissa. They derived a subtle constraint formulation which ensures that the tangent vector to the centerline is continuous and that the length of links remains constant. These two works are worth mentioning since they avoid an empirical tuning of repulsive forces. Yet, [START_REF] Keaveny | Spiral swimming of an artificial micro-swimmer[END_REF] computed the constraint forces and torques with an iterative penalty scheme instead of using an explicit formulation.

Finally, it is worth mentioning that the bead model proposed in [START_REF] Montesi | Brownian dynamics algorithm for bead-rod semiflexible chain with anisotropic friction[END_REF] circumvents the inextensibility difficulty by imposing constraints on the relative velocities of each successive segments, so that their relative distance is kept constant. Using bending potential, [START_REF] Montesi | Brownian dynamics algorithm for bead-rod semiflexible chain with anisotropic friction[END_REF] permit overlap between beads with restoring torque (cf. Fig. 2). A Lagrangian multiplier formulation of tensile forces is also used in [START_REF] Doyle | Dynamic simulation of freely draining flexible polymers in steady linear flows[END_REF], which is equivalent to a prescribed equal distance between successive beads. Again, inextensibility condition does not prevent bead overlapping due to bending in this formulation. The computation of contact forces which is proposed in the following Section 2.2 generalizes the Lagrangian multiplier formulation of [START_REF] Montesi | Brownian dynamics algorithm for bead-rod semiflexible chain with anisotropic friction[END_REF] to generalized forces. Using more complex constraints involving both translational and angular velocities, we show that it is possible to accommodate both non-overlapping and inextensibility conditions without additional repulsive forces (using the rolling no-slip contact with the gears model detailed in Section 2.3). This proposed general formulation is also well suited for any type of kinematic boundary conditions as illustrated in Section 3.4.

Generalized forces, virtual work principle and Lagrange multipliers

The model and formalism proposed in this article rely on earlier work in Analytical Mechanics and Robotics [START_REF] Nikravesh | An overview of several formulations for multibody dynamics[END_REF][START_REF] Joshi | Design and analysis of a spherical mobile robot[END_REF]. The concept of generalized coordinates and constraints which has proven to be very useful in these contexts is described here. Generalized coordinates refer to a set of parameters which uniquely describes the configuration of the system relative to some reference parameters (positions, angles, . . . ). For describing objects of complex shape, let us consider the position r i of each bead i ∈ {1, N b } with associated orientation vector p i which is defined by three Euler angles p ≡ (θ, φ, ψ). In the following, any collection of vector population (r 1 , ..r i , ..r N b ) ≡ R will be capitalized, so that R is a vector in R 3N b . Hence the collection of orientation vectors p i will be denoted P, which is a vector of length 3N b , the collection of velocities dr i dt = ṙi = v i , will be denoted V, the collection of angular velocity ṗi ≡ ω i will be , the collection of forces f i , F, the collection of torques γ i , Γ . All V, , F and Γ are vectors in R 3N b .

Let us then define some generalized coordinate q i for each bead, which is defined by q i ≡ (r i , p i ) ≡ {r 1,i , r 2,i , r 3,i , θ i , φ i , ψ i } so that the collection of generalized positions (q 1 , ..q i , ..q N b ) ≡ Q is a vector in R 6N b . Generalized velocities are then defined by vectors qi ≡ (v i , ω i ) with associated generalized collection of velocities Q.

Articulated systems are generically submitted to constraints which are either holonomic, non-holonomic or both [START_REF] Bailey | Efficient constraint dynamics using MILC SHAKE[END_REF]. Holonomic constraints do not depend on any kinematic parameter (i.e any translational or angular velocity) whereas nonholonomic constraints do.

In the following we consider non-holonomic linear kinematic constraints associated with generalized velocities of the form [START_REF] Greenwood | Classical Dynamics[END_REF] J J J Q + B = 0, [START_REF] Jian | A combined Wormlike-Chain and Bead Model for dynamic simulations of long linear DNA[END_REF] such that J J J is a 6N b × N c matrix and B is a vector of N c components. N c is the number of constraints acting on the N b beads. B and J J J might depend (even non-linearly) on time t and generalized positions Q, but do not depend on any velocity of vector Q, so that relation [START_REF] Jian | A combined Wormlike-Chain and Bead Model for dynamic simulations of long linear DNA[END_REF] is linear in Q. In subsequent sections, we provide specific examples for which this class of constraints are useful. Here we describe, following [START_REF] Nikravesh | An overview of several formulations for multibody dynamics[END_REF][START_REF] Greenwood | Classical Dynamics[END_REF] how such constraints can be handled thanks to some generalized force that can be defined from Lagrange multipliers. The idea formulated to include constraints in the dynamics of articulated systems is to search additional forces which could permit to satisfy these constraints. First, one must rely on generalized forces f i ≡ (f i , γ i ) which include forces and torques acting on each bead, whose collection (f 1 , f i , ..f N b ) is denoted F. Generalized forces are defined such that the total work variation δW is the scalar product between them and the generalized coordinates variations δQ δW = F • δQ = F • δR + Γ • δP, [START_REF] Purcell | Life at low Reynolds number[END_REF] so that, on the right hand side of (6) one also gets the translational and the rotational components of the work. Then, the idea of virtual work principle is to search some virtual displacement δQ that will generate no work, so that

F • δQ = 0. (7)
At the same time, by rewriting [START_REF] Jian | A combined Wormlike-Chain and Bead Model for dynamic simulations of long linear DNA[END_REF] in differential form J J J dQ + Bdt = 0, [START_REF] Fogelson | A fast numerical method for solving the three-dimensional Stokes' equations in the presence of suspended particles[END_REF] admissible virtual displacements, i.e. those satisfying constraints [START_REF] Fogelson | A fast numerical method for solving the three-dimensional Stokes' equations in the presence of suspended particles[END_REF], should satisfy J J J δQ = 0. [START_REF] Stockie | Simulating the motion of flexible pulp fibres using the immersed boundary method[END_REF] Combining the N c constraints ( 9) with ( 7) is possible using any linear combination of these constraints. Such linear combination involves N c parameters, the so-called Lagrange multipliers which are the components of a vector λ in R N c . Then from the difference between [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF] and the N c linear combination of (9) one gets

(F -λ • J J J ) • δQ = 0. ( 10
)
Prescribing an adequate constraint force

F c = λ • J J J , (11)
permits to satisfy the required equality for any virtual displacement. Hence, the constraints can be handled by forcing the dynamics with additional forces, the amplitude of which are given by Lagrange multipliers, yet to be found. Note also, that this first result implies that both translational forces and rotational torques associated with the N c constraints are both associated with the same Lagrange multipliers.

This formalism is particularly suitable for low Reynolds number flows for which translational and angular velocities are linearly related to forces and torques acting on beads by the mobility matrix

M V = M F Γ + V ∞ ∞ + C : E ∞ . ( 12
) V ∞ = (v ∞ 1 , . . . , v ∞ N b ) and ∞ = (ω ∞ 1 , . . . , ω ∞ N b
) correspond to the ambient flow evaluated at the centers of mass r i . E ∞ is the rate of strain 3 × 3 tensor of the ambient flow. C is a third rank tensor called the shear disturbance tensor, it relates the particles velocities and rotations to E ∞ [START_REF] Wajnryb | Generalization of the Rotne-Prager-Yamakawa mobility and shear disturbance tensors[END_REF]. Matrix M (and tensor C) can also be re-organized into a generalized mobility matrix M M M (generalized tensor C C C resp.) in order to define the linear relation between the previously defined generalized velocity and generalized force

Q = M M M F + V V V ∞ + C C C : E ∞ , ( 13
) where V V V ∞ = (v ∞ 1 , ω ∞ 1 , . . . , v ∞ N b , ω ∞ N b ).
The explicit correspondence between the classical matrix M and the hereby proposed generalized coordinate formulation M M M is given in Appendix A. Hence, as opposed to the Euler-Lagrange formalism of classical mechanics, the dynamics of low Reynolds number flows does not involve any inertial contribution, and provide a simple linear relationship between forces and motion. In this framework, it is then easy to handle constraints with generalized forces, because the total force will be the sum of the known hydrodynamic forces F h , elastic forces F e , inner forces associated to active fibers F a and the hereby discussed and yet unknown contact forces F c to verify kinematic constraints

F = F ′ + F c , with (14)
F ′ = F h + F e + F a . ( 15
)
Hence, if one is able to compute the Lagrange multipliers λ, the contact forces will provide the total force by linear superposition [START_REF] Decoene | Microscopic modelling of active bacterial suspensions[END_REF], which gives the generalized velocities with [START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF]. Now, let us show how to compute the Lagrange multiplier vector. Since the generalized force is decomposed into known forces F ′ and unknown contact forces F c = λ • J J J , relations ( 14) and ( 13) can be pooled together yielding

M M M F c = M M M λJ J J = Q -M M M F ′ -V V V ∞ -C C C : E ∞ . ( 16
)
So that, using (5),

J J J M M M J J J T λ = -B -J J J M M M F ′ + V V V ∞ + C C C : E ∞ , (17)
one gets a simple linear system to solve for finding λ, where J J J T stands for the transposition of matrix J J J .

The Gears Model

The Euler-Lagrange formalism can be readily applied to any type of non-holonomic constraint such as (3). In the following, we propose an alternative model based on no-slip condition between the beads: the Gears Model. This constraint, first introduced in a Bead Model (BM) by [START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF], conveniently avoid numerical tricks such as artificial gaps and repulsive forces.

However, [START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF] and [START_REF] Yamamoto | Dynamic simulation of fiber suspensions in shear flow[END_REF] relied on to an iterative procedure to meet requirements. Here, we use the Euler-Lagrange formalism to handle the kinematic constraints associated to the Gears Model.

Considering two adjacent beads (Fig. 4), the velocity v c 1 at the contact point must be the same for each sphere:

v 1 c 1 -v 2 c 1 = 0. (18)
v L c 1 and v R c 1 are respectively the rigid body velocity at the contact point on bead 1 and bead 2. Denote σ 1 the vectorial no-slip constraint. ( 18) becomes

σ 1 (ṙ 1 , ω 1 , ṙ2 , ω 2 ) = 0, (19)
i.e.

[ṙ 1 -

ae 12 × ω 1 ] -[ṙ 2 -ae 21 × ω 2 ] = 0, (20)
where e 12 is the unit vector connecting the center of bead 1, located at r 1 , to the center of bead 2, located at r 2 (e 12 = e 2e 1 ). The orientation p i vector attached to bead i, is not necessary to describe the system. Hence, from (20) one realizes that σ 1 is linear in translational and rotational velocities. Therefore Eq. ( 19) can be reformulated as

σ 1 ( Q) = J 1 Q = 0, (21)
where, Q is the collection vector of generalized velocities of the two-bead assembly

Q = [ṙ 1 , ω 1 , ṙ2 , ω 2 ] T , ( 22
)
J 1 is the Jacobian matrix of σ 1 :

J 1 kl = ∂σ 1 k ∂ Q l , k = 1, . . . , 3, l = 1, . . . , 12, (23)
J 1 = J 1 1 J 1 2 = I 3 -ae × 12 -I 3 ae × 21 , (24)
and

e × =   0 -e 3 e 2 e 3 0 -e 1 -e 2 e 1 0   . ( 25
)
For an assembly of N b beads, N b -1 no-slip vectorial constraints must be satisfied. The Gears Model (GM) total Jacobian matrix J GM is block bi-diagonal and reads

J GM =       J 1 1 J 1 2 J 2 2 J 2 3 . . . . . . J N b -1 N b -1 J N b -1 N b       (26)
where J α β is the 3 × 6 Jacobian matrix of the vectorial constraint α for the bead β.

The kinematic constraints for the whole assembly then read

J GM Q = 0. (27)
The associated generalized forces F c are obtained following Section 2.2.

Elastic forces and torques

We are considering elastohydrodynamics of homogeneous flexible and inextensible fibers. These objects experience bending torques and elastic forces to recover their equilibrium shape. Bending moments derivation and discretization are provided. Then, the role of bending moments and constraint forces is addressed in the force and torque balance for the assembly.

Bending moments

The bending moment of an elastic beam is provided by the constitutive law [START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Bishop | Investigation of bend and shear waves in a geometrically exact elastic rod model[END_REF]

m(s) = K b t × dt ds , (28)
where K b (s) is the bending rigidity, t is the tangent vector along the beam centerline and s is the curvilinear abscissa. Using the Frenet-Serret formula

dt ds = κn, (29)
the bending moment writes

m(s) = K b κb, (30)
where κ(s) defines the local curvature, n(s) and b(s) are the normal and binormal vectors of the Frenet-Serret frame. When the link considered is not straight at rest, with an equilibrium curvature κ eq (s), [START_REF] Schroeder | Effect of hydrodynamic interactions on DNA dynamics in extensional flow: simulation and single molecule experiment[END_REF] is modified into

m(s) = K b κ -κ eq b. (31)
Here, the beam is discretized into N b -1 rigid rods of length l = 2a (cf. Fig. 5). Inextensible rods are made up of two bond beads and linked together by a flexible joint with bending rigidity K b . Bending moments are evaluated at joint locations

s i = (i -1)l for i = 2, . . . , N b -1,
where s i correspond to the curvilinear abscissa of the mass center of the ith bead.

The bending torque on bead i is then given by

γ b i = m(s i+1 ) -m(s i-1 ), (32)
with m(s i ) = K b κ(s i )b(s i ). See Fig. 5 for the torque computation on a beam discretized with four rods. The local curvature κ(s i ) is approximated using the sine rule [START_REF] Lowe | Dynamics of filaments: modelling the dynamics of driven microfilaments[END_REF] κ(s i ) =

1 a 1 + e i-1,i .e i,i+1 2 (33)
where e i-1,i is the unit vector connecting the center of mass of bead i -1 to the center of mass of bead i. This elementary geometric law provides the radius of curvature R(s i ) = 1/κ(s i ) of the circle circumscribing neighboring bead centers r i-1 , r i and r i+1 .

A more general version of the discrete curvature proposed in [START_REF] Fauci | A computational model of aquatic animal locomotion[END_REF] can also be used in the case of three dimensional motion. In that case, the curvature of the fiber is discretized as in [START_REF] Fauci | A computational model of aquatic animal locomotion[END_REF] κ(s i ) =

e i-1,i × e i,i+1 2a , (34)
where, again, e i-1,i is the unit vector connecting the center of mass of bead i -1 to the center of mass of bead i. The bending moment reads m(s i ) = K b κ(s i ). [START_REF] Wada | Non-equilibrium hydrodynamics of a rotating filament[END_REF] To include the effect of torsional twisting about the axis of the fiber, one would have to compute the relative orientation between the frames of reference attached to the beads using Euler angles [START_REF] Keaveny | Spiral swimming of an artificial micro-swimmer[END_REF] (see Section 2.2) or unit quaternions as in [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF]. This would provide the rate of change of the twist angle along the fiber centerline and thus the twisting torque acting on each bead. In the following, only bending effects are considered.

Force and moment acting on each bead

The Gears Model proposed in this paper does not need to consider gaps to allow bending. F c also ensures the connectivity condition and circumvent the use of repulsive forces as distances between adjacent bead surfaces remain constant. More specifically, the tangential components of the force F c , which is only one part of the generalized force F c , acts as tensile force.

For each bead i, contact forces applied from bead i to bead i + 1 at contact point c i between bead i and i + 1 (Fig. 4 for two beads) is denoted f c i . From Newton third law at contact point c i , the contact force applied to bead i from bead i + 1 is obviously -f c i . Total force acting on bead i from contact, and hydrodynamic forces f h i reads

f i = f c i-1 -f c i + f h i (36)
Similarly, the contact force

f c i at point c i produces a moment m c i = at i × f c i associated with local tangent vector t i = e i,i+1
and distance a to the neutral fiber at point c i . Total moment acting on bead i from contact points moments, elastic and hydrodynamic torques are then

γ i = m c i-1 -m c i + γ b i + γ h i . (37)
The contribution of contact force and contact moment acting on bead i exactly equals the contribution of the generalized contact force. Indeed, using the kinematic constraints Jacobian ( 26) in [START_REF] Zhu | Drag of a flexible fiber in a 2D moving viscous fluid[END_REF], and computing the force and torque contributions, one exactly recovers the first and the second contributions of the right-hand-side of ( 36) and [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF]. In Appendix B, we also show that this model is consistent with classical formulation for slender body force and moment balance when the bead radius tends to zero.

Hydrodynamic coupling

Moving objects (rigid or flexible fibers) in a viscous fluid experience hydrodynamic forcing. The interactions are mediated by the fluid flow perturbations which can alter the motion and the deformation of the fibers in a moderately concentrated suspension. The existence of hydrodynamic interactions has also an effect on a single fiber dynamics while different parts of the fiber can respond to the ambient flow but also to local flow perturbations related to the fiber deformation. Resistive Force Theory (RFT) can be used to estimate the fiber response to a given flow assuming that the fiber is modeled by a large series of slender objects [START_REF] Lauga | Floppy swimming: viscous locomotion of actuated elastica[END_REF][START_REF] Coq | Three-dimensional beating of magnetic microrods[END_REF]. Slender body theory has also been used [START_REF] Tornberg | Simulating the dynamics and interactions of flexible fibers in Stokes flows[END_REF][START_REF] Li | The sedimentation of flexible filaments[END_REF] to relate local balance of drag forces with the filament forces upon the fluid resulting in a dynamical system to model the deformation of the fiber centerline. This model provided interesting results on the stretch-coil transition of fibers in vortical flows.

In our beads model, the fiber is composed of spherical particles to account for the finite width of its cross-section. The hydrodynamic interactions are provided through the solution of the mobility problem which relates forces, torques to the translational and rotational velocities of the beads. This many-body problem is non-linear in the instantaneous positions of all particles of the system. Approximate solutions of this complex mathematical problem can be achieved by limiting the mobility matrices to their leading order. The simplest model is called free drain as the mobility matrix is assumed to be diagonal neglecting the HI with neighboring spheres. Pairwise interactions are required to account for anisotropic drag effects within the beads composing the fiber. The Rotne-Prager-Yamakawa (RPY) approximation is one of the most commonly used methods of including hydrodynamic interactions. This widely used approach has been recently updated by Wajnryb et al. [START_REF] Wajnryb | Generalization of the Rotne-Prager-Yamakawa mobility and shear disturbance tensors[END_REF] for the RPY translational and rotational degrees of freedom, as well as for the shear disturbance tensor C which gives the response of the particles to an external shear flow [START_REF] Wagner | The extended finite element method for rigid particles in Stokes flow[END_REF].

Numerical implementation 2.6.1. Integration scheme and algorithm

The kinematics of the constrained system results from the superposition of individual bead motions. Positions are obtained from the temporal integration of the equation of motion with a third order Adams-Bashforth scheme

dr i dt = v i , (38)
where r i , v i are the position and translational velocity of bead i.

The time step t used to integrate [START_REF] Swan | Modeling hydrodynamic self-propulsion with Stokesian Dynamics. Or teaching Stokesian Dynamics to swim[END_REF] is fixed by the characteristic bending time [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF] t < µ(2a) 4 K b , [START_REF] Bilbao | Navigation and chemotaxis of nematodes in bulk and confined fluids[END_REF] where µ is the suspending fluid viscosity.

The evaluation of bead interactions must follow a specific order. Elastic and active forces can be computed in any order.

Constraint forces and torques must be estimated afterwards as they depend on F ′ . Then velocities and rotations are obtained from the mobility relation. And finally, bead positions are updated. [START_REF] Schlagberger | Orientation of elastic rods in homogeneous Stokes flow[END_REF] to get F e , 3. Add active forcing F a and ambient velocity V ∞ if any, 4. Compute the Jacobian matrix associated with non-holonomic constraints J (Q), 5. Solve [START_REF] Simons | The dynamics of sperm detachment from epithelium in a coupled fluid-biochemical model of hyperactivated motility[END_REF] to get the constraint forces F c = λJ , 6. Sum all the forcing terms F = F e + F a + F c , 7. Apply mobility relation [START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF] to obtain the bead velocities Q, 8. Integrate [START_REF] Swan | Modeling hydrodynamic self-propulsion with Stokesian Dynamics. Or teaching Stokesian Dynamics to swim[END_REF] to get the new bead positions.

Implementation of the Joint Model

To provide a comprehensive comparison with previous works, we exploit the flexibility of the Euler-Lagrange formalism to implement the Joint Model as described in [START_REF] Skjetne | Simulation of single fiber dynamics[END_REF] supplemented with hydrodynamic interactions. The joint constraint for two neighboring beads reads ṙi -

(a + ε g )p i × ω i -ṙi+1 + (a + ε g )p i+1 × ω i+1 = 0. (40)
Using the Euler-Lagrange formalism, ( 40) is reformulated with the Joint Model (JM) Jacobian matrix

J JM Q = 0, (41)
where J JM has the same structure as in [START_REF] Gao | Numerical and experimental study of a rotating magnetic particle chain in a viscous fluid[END_REF] and

J i = J i 1 J i 2 = I 3 -(a + ε g )p × i -I 3 -(a + ε g )p × i+1 . (42)
Accordingly, the corresponding set of forces and torques F c are obtained from Section 2.2. As mentioned in Section 2.1, such formulation does not prevent beads from overlapping when bending occurs. A repulsive force F r is added according to [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF] (the force profile proposed by [START_REF] Skjetne | Simulation of single fiber dynamics[END_REF] is very stiff, thus very constraining for the time step):

F r i j =        -F 0 exp(- d i j +δ D d 0 )e i j , d i j ≤ -δ D , -F 0 ( 1 2 - d i j 2δ D )e i j , -δ D < d i j ≤ δ D , 0, r i j > δ D . ( 43
)
δ D is an artificial surface roughness, d i j is the surface to surface distance. d i j < 0 indicates overlapping between beads i and j. d 0 is a numerical damping distance which has to be tuned to prevent overlapping. F 0 is the repulsive force scale chosen in order to avoid numerical instabilities. To deal with this issue, [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF] proposed to evaluate F 0 from bending and viscous stresses. A slight modification of their formula for inertialess particles yields

F 0 = C 1 6π µL v ∞ -v + C 2 K b E b L 3 , (44)
the bar denotes the average over the constitutive beads or joints where C 1 and C 2 are adjustable constants. E b is the bending energy

E b = N b -1 i=1 K b κ(s i ) -κ eq (s i ) 2 . ( 45
)
Bending moments are evaluated at the joint locations Similarly to [START_REF] Schlagberger | Orientation of elastic rods in homogeneous Stokes flow[END_REF], bending torque on bead i is

s J i = (a + ε g ) + (i -1) × 2(a + ε g ), i = 1, . . . , N b -1. Joint curvature is given by κ s J i = 2 a + ε g 1 + p i .p i+1 2 . (46)
γ b i = m s J i -m s J i-1 . (47)
Bead orientation p i is integrated with a third order Adams-Bashforth scheme

dp i dt = ω i × p i . ( 48
)
The procedure is similar to the Gears Model. p i are initialized together with the positions. The repulsive force F r is added to F ′ and can be computed between step 1 and 5 of the aforementioned algorithm. Time integration of ( 48) is performed at step 8.

Constraints and numerical stability

At each time step, the error on kinematic constraints ǫ is evaluated, after application of the mobility relation [START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF], between step 7 and step 8:

ǫ GM (t) = J GM Q 2 = N b -1 i=1 v L c i -v R c i 2 1/2 (49)
for the Gears Model, and

ǫ JM (t) = J JM Q 2 ( 50
)
for the Joint Model.

To verify the robustness of both models and Lagrange formulation, a numerical study is carried out on a stiff configuration.

A fiber of aspect ratio r p = 10 with bending ratio BR = 0.01 is initially aligned with a shear flow of magnitude γ = 5 s -1 . For this aspect ratio, N b = 10 beads are used to model the fiber with the Gears Model.

Joint Model involves additional items to be fixed. N b = 9 spheres are separated by a gap width 2ε g = 0.25a. The repulsive force is activated when the surface to surface distance d i j reaches the artificial surface roughness δ D = 2(a + ε g )/10. The remaining coefficients are set to reduce numerical instabilities without affecting the Physics of the system:

d 0 = (a + ε g )/4
,

C 1 = 5 and C 2 = 0.5.
Fig. 6 shows the evolution of the maximal mean deviation from the no-slip/joint constraint ǭM = max t ǫ(t)/(N b -1)

normalized with the maximal shear velocity γ L depending on the dimensionless time step γ t. First, one can observe that for both Joint and Gears models, ǭM / γ L weakly depends on γ t and the resulting motion of the beads complies very precisely with the set of constraints, within a tolerance close to unit roundoff (< 2.10 -16 ). Secondly, Joint Model is unstable for time steps 100 times smaller than Gears Model. The onset for numerical instability indicates that the repulsive force stiffness dominates over bending, thus dictating and restricting the time step. As a comparison, [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF] matched connectivity constraints within 1% error for each fiber segment. To do so, they had to use an iterative scheme reducing the time step by 1/3 each iteration to meet requirements and limit overlapping between adjacent segments. For similar results, a stiff configuration, such as the sheared fiber, is therefore more efficiently simulated with the Gears Model.

Thirdly, inset of Fig. 6 shows that, for a given time step, the Gears Model constraints ǭM / γ L are satisfied whatever the shear magnitude. Hence, (39) ensures unconditionally numerical stability as bending is the only limiting effect for the Gears Model.

Hence, the robustness of the Euler-Lagrange formalism and the numerical integration we chose provide a strong support to the Gears Model over the Joint Model.

As a final remark to this section, it is important to mention that the numerical cost of the proposed method strongly depends on the choice for the mobility matrix computation, as usual for bead models. If the mobility matrix is computed taking into account full hydrodynamic interactions with Stokesian Dynamics, most of the numerical cost will come from its evaluation in this case. This limitation could be overcame using more sophisticated methods such as Accelerated Stokesian Dynamics [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] or Force Coupling Method [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF]. Moreover, when considering Rotne-Prager-Yamakawa mobility matrix, its cost only requires the evaluation of O ((6N b ) 2 ) terms. Furthermore, the main algorithmic complexity of bead models does not come from the time integration of the bead positions which only requires a matrix-vector multiplication (13) at an O ((6N b ) 2 ) cost. Fast-multipole formulation of a Rotne-Prager-Yamakawa matrix can even provide an O (6N b ) cost for such matrix-vector multiplication [START_REF] Liang | A fast multipole method for the Rotne-Prager-Yamakawa tensor and its applications[END_REF].

The main numerical cost indeed comes from the inversion of the contact forces problem [START_REF] Simons | The dynamics of sperm detachment from epithelium in a coupled fluid-biochemical model of hyperactivated motility[END_REF]. It is worth noting that this linear problem is N c × N c which is slightly different from N b × N b , but of the same order. Furthermore, problem [START_REF] Simons | The dynamics of sperm detachment from epithelium in a coupled fluid-biochemical model of hyperactivated motility[END_REF] gives a direct, single step procedure to compute the contact forces, as opposed to previous other attempts [START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF][START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF][START_REF] Keaveny | Spiral swimming of an artificial micro-swimmer[END_REF] which required iterative procedures to meet forces requirements, involving the mobility matrix inversion at each iteration. The cost for the inversion of ( 17) lies in-between O (N 2 c ) and O (N 3 c ) depending on the inversion method.

Validations

Jeffery orbits of rigid fibers

Much of our current understanding of the behavior of fibers experiencing a shear flow has come from the work of Jeffery [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] who derived the equation for the motion of an ellipsoidal particle in Stokes flow. The same equation can be used for the motion of an axisymmetric particle by using an equivalent ellipsoidal aspect ratio. Rigid fibers can be approximated by elongated prolate ellipsoids. An isolated fiber in simple shear flow rotates in a periodic orbit while the center of mass simply translates in the flow (no migration across streamlines). The period T (51) is a function of the aspect ratio of the fiber and the flow shear rate while the orbit depends on the initial orientation of the object relative to the shear plane T = 2π (r e + 1/r e ) γ .

(

) 51
γ is the shear rate of the carrying flow. r e is the equivalent ellipsoidal aspect ratio which is related to the fiber aspect ratio r p (length of the fiber over diameter of the cross-section which turns out to r p = N b with N b beads). The fiber is initially placed in the plane of shear and is composed on N b beads. No gaps between beads is required in the Joint Model because the fiber is rigid and flexibility deformations are negligible. We have compared the results with two relations for r e : Cox [START_REF] Cox | The motion of long slender bodies in a viscous fluid. Part 2. Shear flow[END_REF] r e = 1.24r p ln(r p ) ,

and Larson [START_REF] Larson | The Structure and Rheology of Complex Fluids[END_REF] r e = 0.7r p .

This classic and simple test case has been successfully validated in [START_REF] Yamamoto | A method for dynamic simulation of rigid and flexible fibers in a flow field[END_REF][START_REF] Yamanoi | Stokesian dynamics simulation of the role of hydrodynamic interactions on the behavior of a single particle suspending in a Newtonian fluid. Part 1. 1D flexible and rigid fibers[END_REF][START_REF] Skjetne | Simulation of single fiber dynamics[END_REF]. Both the Joint and Gears models give a correct prediction of the period of Jeffery orbits (Fig. 7). The scaled γ T of simulations remains within the two evolutions based on Eqs. ( 52) and [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF]. We have tried to compare it with the linear spring model proposed by Gauger and Stark [START_REF] Gauger | Numerical study of a microscopic artificial swimmer[END_REF] (and used by Slowicka et al. [START_REF] Slowicka | Dynamics of fibers in a wide microchannel[END_REF] with a more detailed formulation of hydrodynamic interactions). In this latter model, there is no constraint on the rotation of beads and the simulations failed to reproduce Jeffery orbits (the fiber does not flip over the axis parallel to the flow).

Flexible fiber in a shear flow

The motion of flexible fibers in a shear flow is essential in paper making or composite processing. Prediction and control of fiber orientations and positions are of particular interest in the study of flocks disintegration. Many models have been designed to predict fiber dynamics and much experimental work has been conducted. The wide variety of fiber behaviors depends on the ratio of bending stresses over shear stress, which is quantified by a dimensionless number, the bending ratio BR [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF][START_REF] Forgacs | Particle motions in sheared suspensions: IX. Spin and deformation of threadlike particles[END_REF] BR = E(ln 2r e -1.5)

µ γ 2r 4 p ( 54
)
E is the material Young's Modulus and µ is the suspending fluid viscosity. : theoretical law [START_REF] Hoffman | The dynamics of the coil-stretch transition for long, flexible polymers in planar mixed flows[END_REF] with r e given by ( 53), : theoretical law [START_REF] Hoffman | The dynamics of the coil-stretch transition for long, flexible polymers in planar mixed flows[END_REF] with r e given by ( 52), P: Gears Model, e: Joint Model. In the following, we investigate the response of the Gears Model with known results on flexible fiber dynamics. [START_REF] Forgacs | Particle motions in sheared suspensions: IX. Spin and deformation of threadlike particles[END_REF][START_REF] Forgacs | Particle motions in sheared suspensions: X. Orbits of flexible threadlike particles[END_REF] analyzed the motion of flexible threadlike particles in a shear flow depending on BR. They observed important drifts from the Jeffery orbits and classified them into categories. Yet, the goal of this section is not to carry out an in-depth study on these phenomena. Instead, the objective is to show that our model can reproduce basic features characteristic of sheared flexible filaments. We analyze first the influence of intrinsic deformation on the motion.

Effect of permanent deformation

If a fiber is straight at rest, it will symmetrically deform in a shear flow. When aligned with the compressive axes of the ambient rate of strain E ∞ , the fiber adopts the "S-shape" observed in Fig. 8(a). When aligned with the extensional axes, tensile forces turn the rod back to its equilibrium shape. This symmetry is broken when the filament is initially slightly deformed or has a permanent deformation at rest, i.e. a nonzero equilibrium curvature κ eq > 0. An initial small perturbation of the shape of a straight filament aligned with flow can induce large deformations during the orbit. This phenomenon is known as the buckling instability whose onset and growth are quantified with BR [START_REF] Becker | Instability of elastic filaments in shear flow yields first-normal-stress differences[END_REF][START_REF] Guglielmini | Buckling transitions of an elastic filament in a viscous stagnation point flow[END_REF]. Fig. 8(b) illustrates the evolution of a flexible sheared filament with BR = 0.04 and a very small intrinsic deformation κ eq = 1/(100L). The equilibrium dimensionless radius of curvature is 2R eq /L = 200. During the tumbling motion it decreases to a minimal value of 2R min /L = 0.26. Buckling thus increases by 770 times the maximal fiber curvature.

Maximal fiber curvature

[74] measured the radius of curvature R of sheared fiber for aspect ratios r p ranging from 283 to 680. They reported on the evolution of the minimal value R min , i.e. the maximal curvature κ max , with BR. [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF] used the Joint Model with prolate spheroids but no hydrodynamic interactions and compared their results with [START_REF] Salinas | Bending and breaking fibers in sheared suspensions[END_REF]. Both experimental results from [START_REF] Salinas | Bending and breaking fibers in sheared suspensions[END_REF] and simulations from [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF] are accurately reproduced by the Gears Model. Hydrodynamic interactions between fiber elements play an important role in the bending of flexible filaments [START_REF] Lindström | Simulation of the motion of flexible fibers in viscous fluid flow[END_REF][START_REF] Slowicka | Dynamics of fibers in a wide microchannel[END_REF][START_REF] Salinas | Bending and breaking fibers in sheared suspensions[END_REF]. As mentioned in Section 2.5 the use of spheres to build any arbitrary object is well suited to compute these hydrodynamic interactions. However, modeling rigid slender bodies in a strong shear flow becomes costly when increasing the fiber aspect ratio. First, the aspect ratio of a fiber made up of N b spheres is r p = N b . Each time iteration requires the computation of !: current simulations with aspect ratio r p = 35 and intrinsic curvature κ eq = 0; ": current simulations with aspect ratio r p = 35 and intrinsic curvature κ eq = 1/(10L); simulation results from [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF] with κ eq = 1/(10L): (E: r p = 50, P: r p = 100, e: r p = 150, 1: r p = 280); +: experimental measurements from [START_REF] Salinas | Bending and breaking fibers in sheared suspensions[END_REF], r p = 283.

M and C : E ∞ and the inversion of a linear system (17) corresponding to N c relations of constraints with

N c ≥ 3(N b -1).
Secondly, for a given shear rate γ and bending ratio BR, Young's modulus increases as r 4

p . According to [START_REF] Bilbao | Navigation and chemotaxis of nematodes in bulk and confined fluids[END_REF], the time step becomes very small for large E. [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF] partially avoided this issue by neglecting pairwise hydrodynamic interactions (M is diagonal), and by assembling prolate spheroids of aspect ratios r e ∼ 10.

Yet, it is shown in Fig. 9(a) that for a fixed BR, 2R min /L converges asymptotically to a constant value with r p . An asymptotic regime (relative variation less than 2%) is reached for r p ≥ 25. Choosing r p = 35 thus enables a valid comparison with previous results.

Our simulation results compare well with the literature data (Fig. 9(b)) and better match with to experiments than [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF].

When BR ≥ 0.04, the Gears Model clearly underestimates measurements for κ eq = 1/(10L) and overestimates them for κ eq = 0. However, Salinas and Pittman [74] indicated that the error quantification on parameters and measurements is diffi- cult to estimate as the fibers were hand-drawn. Notably, drawing accuracy decreases for large radii of curvature, which leads to the conclusion that the hereby observed discrepancy might not be critical. They did not report the value of permanent deformation κ eq for the fibers they designed, whereas, as evidenced by [71], it has a strong impact on R min . A numerical study of this dependence should be conducted to compare with [71, Fig. 7].

[46] used the same approach as [START_REF] Schmid | Simulations of fiber flocculation: effects of fiber properties and interfiber friction[END_REF] with hydrodynamic interactions to repeat numerically the experiments from [START_REF] Salinas | Bending and breaking fibers in sheared suspensions[END_REF]; but their results, though reliable, were displayed such that direct comparison with previous work is not possible.

To conclude, it should be noted that, in [53], the aspect ratio does not affect 2R min /L for a fixed BR, confirming the asymptotic behavior observed in Fig. 9(a).

Settling fiber

Consider a fiber settling under constant gravity force F ⊥ = F ⊥ e ⊥ acting perpendicularly to its major axis. The dynamics of the system depends on three competing effects: the elastic stresses which tend to return the object to its equilibrium shape, the gravitational acceleration which uniformly translates the object and the hydrodynamic interactions which creates local drag along the filament. After a transient regime, the filament reaches steady state and settles at a constant velocity with a fixed shape (see Figs. 10(a) and 10(b)). This steady state depends on the elasto-gravitational number [START_REF] Keaveny | Dynamics of structures in active suspensions of paramagnetic particles and applications to artificial micro-swimmers[END_REF][START_REF] Cosentino Lagomarsino | Hydrodynamic induced deformation and orientation of a microscopic elastic filament[END_REF] and [START_REF] Li | The sedimentation of flexible filaments[END_REF] examined the contribution of each competing effect by measuring the normal deflection A, i.e. the distance between the uppermost and the lowermost point of the filament along the direction of the applied force (Fig. 10(b)); and the normal friction coefficient γ ⊥ /γ 0 ⊥ as a function of B. γ 0 ⊥ is the normal friction coefficient of a rigid rod. To compute hydrodynamic interactions [START_REF] Cosentino Lagomarsino | Hydrodynamic induced deformation and orientation of a microscopic elastic filament[END_REF] used Stokeslet; [START_REF] Keaveny | Dynamics of structures in active suspensions of paramagnetic particles and applications to artificial micro-swimmers[END_REF], the Force Coupling Method (FCM) [START_REF] Yeo | Simulation of concentrated suspensions using the force-coupling method[END_REF]; [START_REF] Li | The sedimentation of flexible filaments[END_REF], Slender Body Theory.

B = F ⊥ L/K b . (55)
Similar simulations were carried out with both the Joint Model described in Section 2.6.2 and the Gears Model. Fiber of length L = 68a is made out of N b = 31 beads with gap width 2ε g = 0.2a for the Joint Model and N b = 34 for the Gears Model. To avoid both overlapping and numerical instabilities with the Joint Model, the following repulsive force coefficients were selected: Fig. 11 shows that our simulations agree remarkably well with previous results except slight differences with [START_REF] Li | The sedimentation of flexible filaments[END_REF] in the linear regime B < 100. Using Slender Body Theory, [START_REF] Li | The sedimentation of flexible filaments[END_REF] made the assumption of a spheroidal filament instead of a cylindrical one, with aspect ratio r p = 100, i.e. three times larger than other simulations, whence such discrepancies. The Fig. 10. Shape of settling fiber for B = 10 000 in the frame moving with the center of mass (x c , z c ). (a) Metastable "W" shape, t = 12L/V s . (b) Steady "horseshoe" shape at t = 53L/V s . V s is the terminal settling velocity once steady state is reached. normal coefficient (Fig. 11(b)), resulting from hydrodynamic interactions, perfectly matches the value obtained by [START_REF] Keaveny | Dynamics of structures in active suspensions of paramagnetic particles and applications to artificial micro-swimmers[END_REF] with the Force Coupling Method. The FCM is known to better describe multibody hydrodynamic interactions. Such a result thus supports the use of the simple Rotne-Prager-Yamakawa tensor for this hydrodynamic system.

d 0 = (a + ε g )/4, δ D = (a + ε g )/5, C 1 = 0.
Differences between Gears and Joint Models implemented here are quantified by measuring the relative discrepancies on the vertical deflection

A ǫ A = A G -A J A G , (56)
and on the normal friction coefficient γ ⊥ /γ 0

⊥ ǫ γ ⊥ = (γ ⊥ /γ ⊥ 0 ) G -(γ ⊥ /γ ⊥ 0 ) J (γ ⊥ /γ ⊥ 0 ) G . (57)
Discrepancies between Joint and Gears models remain below 5% except at the threshold of the non-linear regime (B ≈ 100) where ǫ A reaches 15% and ǫ γ ⊥ ≈ 7.5%.

In accordance with [START_REF] Cosentino Lagomarsino | Hydrodynamic induced deformation and orientation of a microscopic elastic filament[END_REF], a metastable "W" shape is reached for B > 3000 (Fig. 10(a)) until it converges to the stable "horseshoe" state (Fig. 10(b)).

Actuated filament

The goal of the following sections is to show that the model we proposed is not only valid for passive objects but also for active ones. Elastohydrodynamics also concern swimming at low Reynolds number [START_REF] Purcell | Life at low Reynolds number[END_REF]. Many type of micro-swimmers have been studied both from the experimental and theoretical point of view. Among them two categories are distinguished: ciliates and flagellates. Ciliates propel themselves by beating arrays of short hairs (cilia) on their surface in a synchronized way (Opalina, Volvox, Paramecia). Flagellates undulate and/or rotate their external appendage to push (pull) the fluid from their aft (fore) part (spermatozoa, Chlamydomonas Rheinardtii, Bacillus Subtilis, Eschericia Coli). Recent advances in nanotechnologies allows researchers to design artificial swimming micro-devices inspired by low Reynolds number fauna [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF][START_REF] Zhang | Characterizing the swimming properties of artificial bacterial flagella[END_REF][START_REF] Keaveny | Optimization of chiral structures for microscale propulsion[END_REF].

In that scope, the study of bending wave propagation along passive elastic filament has been investigated by [START_REF] Wiggins | Trapping and wiggling: elastohydrodynamics of driven microfilaments[END_REF][START_REF] Yu | Experimental investigations of elastic tail propulsion at low Reynolds number[END_REF] and [START_REF] Coq | Helical beating of an actuated elastic filament[END_REF][START_REF] Coq | Rotational dynamics of a soft filament: wrapping transition and propulsive forces[END_REF].

The experiment of [START_REF] Yu | Experimental investigations of elastic tail propulsion at low Reynolds number[END_REF] consists in a flexible filament tethered and actuated at its base. The base angle was oscillated sinusoidally in plane with an amplitude α 0 = 0.435 rad and frequency ζ .

Deformations along the tail result from the competing effects of bending and drag forces acting on it. A dimensionless quantity called the Sperm number compares the contribution of viscous stresses to elastic response [START_REF] Wiggins | Flexive and propulsive dynamics of elastica at low Reynolds number[END_REF]

Sp = L ζ (γ ⊥ /L) K b 1/4 = L l ζ . ( 58
)
γ ⊥ is the normal friction coefficient. When using Resistive Force Theory, (γ ⊥ /L) is changed into a drag per unit length coefficient ξ ⊥ . l ζ can be seen as the length scale at which bending occurs. Sp 1 corresponds to a regime at which bending dominates over viscous friction: the whole filament oscillates rigidly in a reversible and symmetrical way. Sp ≫ 1 corresponds to a regime at which bending waves are immediately damped and the free end is motionless [START_REF] Wiggins | Flexive and propulsive dynamics of elastica at low Reynolds number[END_REF]. The experiment of [START_REF] Coq | Rotational dynamics of a soft filament: wrapping transition and propulsive forces[END_REF] is similar to [START_REF] Yu | Experimental investigations of elastic tail propulsion at low Reynolds number[END_REF] except for that the actuation at the base is rotational. Here, the filament was rotated at a frequency ζ forming a base angle α 0 = 0.262 rad with the rotation axis.

In both contributions, the resulting fiber deformations were measured and compared to Resistive Force Theory for several values of Sp. Simulations of such experiments [START_REF] Yu | Experimental investigations of elastic tail propulsion at low Reynolds number[END_REF][START_REF] Coq | Rotational dynamics of a soft filament: wrapping transition and propulsive forces[END_REF] were performed with the Gears Model.

Numerical setup and boundary conditions at the tethered base element

Corresponding kinematic boundary conditions for BM are prescribed with the constraint formulation of the Euler-Lagrange formalism.

Planar actuation

In the case of planar actuation [START_REF] Yu | Experimental investigations of elastic tail propulsion at low Reynolds number[END_REF], we consider that the tethered, i.e. the first, bead is placed at the origin and has no degree of freedom

ṙc 1 = 0, ω c 1 = 0.

Helical actuation

In the case of helical beating [START_REF] Coq | Helical beating of an actuated elastic filament[END_REF][START_REF] Coq | Rotational dynamics of a soft filament: wrapping transition and propulsive forces[END_REF], the anchor point of the filament is slightly off-centered with respect to the rotation axis e x [START_REF] Coq | Helical beating of an actuated elastic filament[END_REF]: r(0) = δ 0 (cf. Fig. 13, left inset). [START_REF] Coq | Helical beating of an actuated elastic filament[END_REF] measured a value δ 0 = 2 mm with a filament length varying from L = 2 cm to 10 cm. Here we take δ 0 = δ0 sin α 0 with δ0 = 2.7a and vary the filament length by changing the number of beads N b to match the experimental range δ 0 /L = 0.1 → 0.02. The position of bead 1 must then follow r c 1 (t) =    δ0 cos(α 0 sin(ζ t)) cos(α 0 cos(ζ t)) δ0 cos(α 0 sin(ζ t)) sin(α 0 cos(ζ t)) δ0 sin(α 0 sin(ζ t))    .

(

) 62
The translational velocity of the first bead ṙ1 (t) is thus constrained by taking the derivative of ( 62)

ṙc 1 (t) =        δ0 α 0 ζ [-cos(ζ t) sin(α 0 sin(ζ t)) cos(α 0 cos(ζ t)) + sin(ζ t) sin(α 0 cos(ζ t)) cos(α 0 sin(ζ t))] δ0 α 0 ζ [-cos(ζ t) sin(α 0 sin(ζ t)) sin(α 0 cos(ζ t)) -sin(ζ t) cos(α 0 cos(ζ t)) cos(α 0 sin(ζ t))] δ0 α 0 ζ cos(ζ t) cos(α 0 sin(ζ t))        . ( 63
)
And the rotational velocity is set to zero ω 1 = 0.

The velocity of the second bead ṙc 2 (t) is prescribed in synchrony with bead 1:

ṙc 2 (t) =        ( δ0 + 2a)α 0 ζ [-cos(ζ t) sin(α 0 sin(ζ t)) cos(α 0 cos(ζ t)) + sin(ζ t) sin(α 0 cos(ζ t)) cos(α 0 sin(ζ t))] ( δ0 + 2a)α 0 ζ [-cos(ζ t) sin(α 0 sin(ζ t)) sin(α 0 cos(ζ t)) -sin(ζ t) cos(α 0 cos(ζ t)) cos(α 0 sin(ζ t))] ( δ0 + 2a)α 0 ζ cos(ζ t) cos(α 0 sin(ζ t))        . ( 64
)
The rotational velocity ω 2 is consistently constrained by the no-slip condition. The three-dimensional curvature κ is discretized with [START_REF] Yamanoi | Stokesian dynamics simulation of the role of hydrodynamic interactions on the behavior of a single particle suspending in a Newtonian fluid. Part 1. 1D flexible and rigid fibers[END_REF].

In both cases, imposing actuation at the base of the filament therefore requires the addition of three vectorial kinematic constraints, [START_REF] Joshi | Design and analysis of a spherical mobile robot[END_REF] and [START_REF] Yamamoto | Dynamic simulation of fiber suspensions in shear flow[END_REF], to the no-slip conditions: N c = 3(N b -1) + 3 × 3. The additional Jacobian matrix J act writes

J act =   I 3 0 3 0 3 0 3 • • • 0 3 0 3 0 3 I 3 0 3 0 3 • • • 0 3 0 3 0 3 0 3 I 3 0 3 • • • 0 3 0 3   . ( 65
)
The corresponding right-hand side B act contains the imposed velocities

B act =   0 0 -ṙ c 2   (66)
for planar beating, and

B act =   -ṙ c 1 0 -ṙ c 2   (67)
for helical beating.

J act and B act are simply appended to J J J and B respectively; corresponding forces and torques F c are computed as explained before in Section 2.2.

Comparison with experiments and theory

The dynamics of the system can be described by balancing elastic stresses (flexion and tension) with viscous drag. Subsequent coupled non-linear equations can be linearized with the approximation of small deflections or solved with an adaptive integration scheme [START_REF] Lauga | Floppy swimming: viscous locomotion of actuated elastica[END_REF][START_REF] Gadêlha | Nonlinear instability in flagellar dynamics: a novel modulation mechanism in sperm migration?[END_REF][START_REF] Camalet | Generic aspects of axonemal beating[END_REF]. [START_REF] Yu | Experimental investigations of elastic tail propulsion at low Reynolds number[END_REF] considered both linear and non-linear theories and included the effect of a sidewall by using the corrected RFT coefficients of [START_REF] Mestre | Low-Reynolds-number translation of a slender cylinder near a plane wall[END_REF].

Planar actuation

Simulations are in good agreement with experiments, linear and non-linear theories for Sp = 1.73, 2.2, and 3.11 (Fig. 12). Even though sidewall effects were neglected here, the Gears Model provides a good description of non-linear dynamics of an actuated filament in Stokes flow. E: experiment, ": Gears Model with no anchoring distance δ 0 /L = 0, 2: Gears Model with δ 0 /L = 0.1 → 0.02 as in [START_REF] Coq | Helical beating of an actuated elastic filament[END_REF].

Helical actuation

Once steady state was reached, [START_REF] Coq | Rotational dynamics of a soft filament: wrapping transition and propulsive forces[END_REF] measured the distance of the tip of the rotated filament to the rotation axis d = r(L) (cf. Fig. 13, left inset). Fig. 13 compares their measures with our numerical results. Insets show the evolution of the filament shape with Sp. The agreement is quite good. Numerical simulations slightly overestimate d for 30 < Sp 4 < 90. This may be due to the lack of information to reproduce experimental conditions and/or to measurement errors. As stated in [START_REF] Coq | Helical beating of an actuated elastic filament[END_REF], taking the anchoring distance δ 0 into account is important to match the low Sperm number configurations where δ 0 /L is non-negligible and the filament is stiff. If the anchoring point was aligned with the rotation axis (δ 0 = 0), the distance to the axis of the rod free end would be d/L = sin α 0 = 0.259 for small Sp, as shown on Fig. 13.

Planar swimming nematode

Locomotion of the nematode Caenorhabditis Elegans is addressed here as its dynamics and modeling are well documented [START_REF] Bilbao | Navigation and chemotaxis of nematodes in bulk and confined fluids[END_REF][START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF]. C. Elegans swims by propagating a contraction wave along its body length, from the fore to the aft (Fig. 14(a)).

Modeling such an active filament in the framework of BM just requires the addition of an oscillating driving torque γ D (s, t) to mimic the internal muscular contractions. To do so, [START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF] used the preferred curvature model. In this model, the driving torque results from a deviation in the centerline curvature from κ D (s, t) = -κ D 0 (s) sin(ks -2π f t), [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] where κ D 0 (s) is prescribed to reproduce higher curvature near the head:

κ D 0 (s) = K 0 , s ≤ 0.5L 2K 0 (L -s)/L, s > 0.5L. ( 69
)
The amplitude K 0 , wave number k and the associated Sperm number Sp = L f µ/K b 1/4 [START_REF] Forgacs | Particle motions in sheared suspensions: IX. Spin and deformation of threadlike particles[END_REF] were tuned to reproduce the measured curvature wave of the free-swimming nematode. They obtained the following set of numerical values: K 0 = 8.25/L, k = 1.5π /L and Sp * = 22.6 1/4 . The quantity of interest to compare with experiments is the distance the nematode travels per stroke V /( f L). K b is assumed to be constant along s and is deduced from the other parameters. As for [START_REF] Schlagberger | Orientation of elastic rods in homogeneous Stokes flow[END_REF], the torque applied on bead i results from the difference in active bending moments across neighboring links

γ D i (t) = m D (s i+1 , t) -m D (s i-1 , t), ( 71
) with m D (s i , t) = K b κ D (s i , t)b(s i ).
γ D is added to F a at step 3 of the algorithm in Section 2.6.1.

To match the aspect ratio of C. Elegans, r p = 16, [START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF] put N b = 15 beads together separated by gaps of width 2ǫ g = 0.2a. Here we assemble N b = 16 beads, avoiding the use of gaps, and employ the same target-curvature wave and numerical coefficient values.

The net translational velocity V * = V /( f L) = 0.0662 obtained with our model matches remarkably well with the numerical results V /( f L) = 0.0664 [START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF] and experimental measurements V /( f L) ≈ 0.07 [START_REF] Bilbao | Navigation and chemotaxis of nematodes in bulk and confined fluids[END_REF].

Cooperative swimming

One of the configurations explored in [START_REF] Llopis | Cooperative motion of intrinsic and actuated semiflexible swimmers[END_REF] has been chosen as a test case for the interactions between in-phase or out-of phase swimmers. Two identical, coplanar C. Elegans swim in the same direction with a phase difference φ which is introduced in the target curvature, and thus in the driving torque, of the second swimmer κ D,2 (s, t) = -κ D 0 (s) sin(ks -2π f t + φ). The initial shape of the swimmers is taken from their steady state. We define d as the distance between their center of mass at initial time (see inset of Fig. 14(b)). Similarly to [83, Fig. 3], our results (Fig. 14(b)) show that antiphase beating enhance the propulsion, whereas in-phase swimming slows the system as swimmers get closer. Even though the model swimmer here is different, the quantitative agreement with [START_REF] Llopis | Cooperative motion of intrinsic and actuated semiflexible swimmers[END_REF] is strikingly good. Numerical work by [START_REF] Fauci | Interaction of oscillating filaments: a computational study[END_REF] also revealed that the average swimming speed of infinite sheets in finite Reynolds number flow is maximized when they beat in opposite phase. The conclusion that closer swimmers do not necessarily swim faster than individual ones has also been reported in [START_REF] Immler | By hook or by crook? Morphometry, competition and cooperation in rodent sperm[END_REF]. They measured a decrease in the swimming speed of 25% for groups of house mouse sperm, as obtained on Fig. 14(b) for d/L = 0.2.

Spiral swimming

Many of the flagellate micro-organisms such as spermatozoa, bacteria or artificial micro-devices use spiral swimming to propel through viscous fluid. Propulsion with rotating rigid or flexible filaments has been thoroughly investigated in the past years [START_REF] Coq | Helical beating of an actuated elastic filament[END_REF][START_REF] Manghi | Propulsion with a rotating elastic nanorod[END_REF][START_REF] Keaveny | Spiral swimming of an artificial micro-swimmer[END_REF][START_REF] Keaveny | Optimization of chiral structures for microscale propulsion[END_REF][START_REF] Cortez | The method of regularized Stokeslets in three dimensions: analysis, validation, and application to helical swimming[END_REF][START_REF] Qian | Shape transition and propulsive force of an elastic rod rotating in a viscous fluid[END_REF][START_REF] Hsu | A 3d motile rod-shaped monotrichous bacterial model[END_REF]. In this section we illustrate the versatility of the proposed model by investigating the effect of the Sperm number and the eccentricity of the swimming gait on the swimming speed of C. Elegans.

Numerical configuration

The target curvature of C. Elegans κ D remains unchanged except for that it is now directed along two components which are orthogonal to the helix axis. A phase difference φ = π/2 is introduced between these two components. The resulting driving moment writes: m D (s i , t) = αK b κ D (s i , t)e ⊥ + β K b κ D (s i , t, φ = π/2)e b . [START_REF] Guglielmini | Buckling transitions of an elastic filament in a viscous stagnation point flow[END_REF] {e , e ⊥ , e b } are body fixed orthonormal vectors. e is directed along the axis of the helix, e ⊥ is a perpendicular vector and e b is the binormal vector completing the basis (Fig. 15 inset). The magnitude of the curvature wave along e ⊥ (resp. e b ) is weighted by a coefficient α (resp. β). The trajectory of a body element in the plane {e ⊥ , e b } describes an ellipse whose eccentricity depends on the value of the ratio β/α. When β/α = 0 the driving torque is two-dimensional and identical to the one used in Section 3.5. When β/α = 1 the magnitude of the driving torque is equal in both direction, the swimming gait describes a circle in the plane {e ⊥ , e b } (see Fig. 15 inset). For the sake of simplicity, here we take {e , e ⊥ , e b } = {e x , e y , e z }.

As in Section 3.4, the curvature is evaluated with [START_REF] Yamanoi | Stokesian dynamics simulation of the role of hydrodynamic interactions on the behavior of a single particle suspending in a Newtonian fluid. Part 1. 1D flexible and rigid fibers[END_REF]. In the following, α = 1 and only β is varied in the range [0; 1].

Results

Fig. 15 compares the planar swimming speed of C. Elegans V * with its "helical" version V depending on the Sperm number defined in Section 3.5 [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] and on the ratio β/α. The Sperm number Sp lies in the range [17 reaches a value such that the swimmer experiences a change in shape which is not helical. This sudden change in shape breaks any periodical motion and makes irrelevant the measurement of a net translational motion. Such limitation is only linked to the choice of the numerical coefficients of the target curvature model.

For the characteristic value Sp * = 22.6 1/4 chosen by [START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF], the purely helical motion provides a swimming speed four times faster than planar beating. Even though the model swimmer is different here, this result qualitatively agrees with the observation of [START_REF] Keaveny | Spiral swimming of an artificial micro-swimmer[END_REF] for which spiral swimming was faster than planar beating. Beyond a critical value Sp ≈ 2.6, planar beating is faster. For β/α = 0.5 the swimming speed is always smaller than for planar beating except when Sp < 2.15.

This last observation is not intuitive. A more extensive study on the effect of the eccentricity of the swimming gait on the swimming speed would be of interest.

Conclusions

We have provided a simple general theoretical framework for kinematic constraints to be used in three-dimensional BM. This framework permits to handle versatile and complex kinematic constraints between flexible assembly of spheres, and/or more complex non-deformable objects at low Reynolds numbers. Using Stokes linearity, this formulation requires, at each time step, the inversion of an O (N c × N c ) linear system for an assembly having N c constraints. Constraints are exactly matched (up to machine precision) and their evaluation is insensitive to time-step. Furthermore, since the formulation explicitly handles mobility matrices, it can be used with any approximation for hydrodynamics interactions, from free drain (no HI) to full Stokesian Dynamics. The proposed framework also implicitly incorporates the generic influence of external flows on kinematic constraints, as opposed to previous BM formulation which necessitates some adjustments to the ambient flow in most cases.

We also propose a simple Gears Model to describe flexible objects, and we show that such model successfully predicts the fiber dynamics in an external flow, its response to an external mechanical forcing and the motion of internally driven swimmers. Quantitative agreement with previous works is obtained for both slender objects (fibers, actuated filaments) and non-slender swimmers (C. Elegans), allowing its use in a wide variety of contexts. The Gears Model is easy to implement and it fulfills several important improvements over previous BM:

• There is no limitation on the fiber curvature, since Gears Model does not need any repulsive force nor gap width to be defined.

• Gears Model is more generic than previous ones, since there is no need for numerical parameter to be tuned.

• When compared with Lagrange multiplier formulation of Joint Model, Gears Model is also much more stable by two orders of magnitude in time-step, a drastic improvement which offers nice prospects for the modeling of complex flexible assemblies.

Finally it should be noted that even if we only consider simple collections of spheres, any complex assembly can be easily treated within a similar framework, which provide interesting prospects in the future modeling of complex micro-organisms, membranes or cytoskeleton micro-mechanics.

Fig. 1 .

 1 Fig. 1. Spring Model: linear spring to keep constant the inter-particle distance.

Fig. 2 .

 2 Fig. 2. Joint Model: overlapping due to bending if no gap between beads.

Fig. 3 .

 3 Fig. 3. Joint Model: c 1 is separated by a gap ε g from the beads.

Fig. 4 .

 4 Fig. 4. Gears Model: contact velocity must be the same for each bead (no-slip condition).

Fig. 5 .

 5 Fig. 5. Beam discretization and bending torques computation of beads 1, 3 and 5. Remaining torques are accordingly obtained: γ b 2 = m(s 3 ) and γ b 4 = -m(s 3 ).

Fig. 6 .

 6 Fig.6. Dependence of the constraints ǭM / γ L on the time step γ t, +: Gears Model, 1: Joint Model. Inset: ǭM / γ L with the Gears Model for a fixed time step given by[START_REF] Bilbao | Navigation and chemotaxis of nematodes in bulk and confined fluids[END_REF] for different values of γ .

Fig. 7 .

 7 Fig.7. Tumbling period T depending on fiber aspect ratio r p .

Fig. 8 .

 8 Fig. 8. Orbit of a flexible filament in a shear flow with BR = 0.04. Temporal evolution is shown in the plane of shear flow. (a) Symmetric "S-shape" of a straight filament, κ eq = 0. (b) Buckling of a permanently deformed rod with an intrinsic curvature κ eq = 1/(100L).

Fig. 9 .

 9 Fig. 9. (a) Minimal radius of curvature depending on fiber length for bending ratios.

  01 and C 2 = 0.01. No adjustable parameters are required for Gears Model.

Fig. 11 .

 11 Fig. 11. (a) Scaled vertical deflection A/L depending on B.: Gears Model, : Joint Model, : FCM results from[START_REF] Keaveny | Dynamics of structures in active suspensions of paramagnetic particles and applications to artificial micro-swimmers[END_REF], : Stokeslets results from[START_REF] Cosentino Lagomarsino | Hydrodynamic induced deformation and orientation of a microscopic elastic filament[END_REF],: Slender body theory results from[START_REF] Li | The sedimentation of flexible filaments[END_REF]. (b) Normal friction coefficient vs. B. : Gears Model, : Joint Model, : FCM results from[START_REF] Keaveny | Dynamics of structures in active suspensions of paramagnetic particles and applications to artificial micro-swimmers[END_REF], : Stokeslets results from[START_REF] Cosentino Lagomarsino | Hydrodynamic induced deformation and orientation of a microscopic elastic filament[END_REF].

( 59 )

 59 Denote α 0 the angle formed between e x and e 1,2 .The trajectory of bead 2 must follow

( 60 )

 60 The translational velocity of the second bead ṙ2 (t) is thus constrained by taking the derivative of (60)ṙc 2 (t) =   -2aα 0 ζ cos(ζ t) sin(α 0 sin(ζ t)) 0 2aα 0 ζ cos(ζ t) cos(α 0 sin(ζ t))   .

Fig. 12 .

 12 Fig. 12. Comparison with experiments and numerical results from [79]. Gears Model results are superimposed on the original Fig. 3 of [79]. Snapshots are shown for four equally spaced intervals during the cycle for one tail with α 0 = 0.435 rad. : experiment,

Fig. 13 .

 13 Fig. 13. Comparison with experiments from [80]. (Insets) Evolution of the filament shape with Sp 4 . Snapshots are shown for twenty equally spaced intervals during one period at steady state. Gray level fades as time progresses. Left inset: δ 0 /L is the distance of the tethered bead to the rotation axis, d/L is the distance of the free end to the rotation axis. (Main figure) Distance of the rod free end to the rotation axis normalized by the filament length d/L.

Fig. 14 .

 14 Fig. 14. (a) Simulated wave motion of a swimming model C. Elegans. The nematode swims leftward and gray level fades as time progresses. Motion is shown in a frame moving with the micro-swimmer center of mass. (b) (Inset) Two C. Elegans beating in the same plane at a distance d in opposite phase ( φ = π ). Nematodes swim leftward and gray level fades as time progresses. (Main figure) Swimming speed of the center of mass of the system V normalized by the isolated swimming speed of C. Elegans V * . !: in-phase motion ( φ = 0); 1: antiphase motion ( φ = π ).

Fig. 15 .

 15 Fig. 15. Helical swimming of C. Elegans. (Inset) Snapshot for Sp * = 22.6 1/4 and β/α = 1.

•

  Initialization: positions r i (0), • Time Loop 1. Evaluate mobility matrix M (Q) and C : E ∞ (see Section 2.5), 2. Calculate local curvatures (33) and bending torques γ b i

  1/4 ; 10001/4 ] = [2.03; 5.62]. The lower bound is dictated by the stability of the helical swimming. When Sp < 2.03, the imposed curvature

Acknowledgements

The authors would like to acknowledge support from the ANR project MOTIMO: Seminal Motility Imaging and Modeling (grant no. ANR-11-MONU-009-01). We are also thankful to Pr Pierre Degond, Yuan-nan Young and Dr Eric E. Keaveny for fruitful discussions.

Appendix A. Correspondence between M and M M M

Matrix M M M defined in (13) results from the rearrangement of the well-known mobility matrix M. This operation is necessary in order to combine constraint Eq. ( 5) and mobility relation [START_REF] Wagner | The extended finite element method for rigid particles in Stokes flow[END_REF] to obtain the constraint forces F c .

Matrix M relates the collection of velocities V = (v 1 , . . . , v N b ) and rotations = (ω 1 , . . . , ω N b ) to the collection of forces F = (f 1 , . . . , f N b ) and torques Γ = (γ 1 , . . . , γ N b )

where M VF is the 3N b × 3N b matrix relating all the bead velocities to the forces applied to their center of mass

to obtain a mobility equation suited for the Euler-Lagrange formalism

) is strictly equivalent to (13).

Appendix B. Asymptotic limit of force and moment balance on the Gears Model

In this appendix we show that slender body formulation for elastic fibers, when applied to Gears Model is consistent with the discrete formulation of force and moments balance [START_REF] Manghi | Propulsion with a rotating elastic nanorod[END_REF] and (37) in the asymptotic limit of small beads.

The force balance equation for a beam is [START_REF] Landau | Theory of Elasticity[END_REF] ∂n is ∂s

n is (s) is the resultant internal stress on a cross-section S(s) at arclength position s along the centerline

for which the tangent vector to neutral fiber centerline is also the unit normal vector to cross-section S(s). f is the force per unit length which contains any supplementary contribution to the internal elastic response of the material (e.g. hydrodynamic force per unit length). The moment balance reads [START_REF] Landau | Theory of Elasticity[END_REF] ∂m is ∂s

where m is (s) is the moment of the flexion and torsion stresses on the cross-section which are related to the local deformation of Frenet-Serret coordinates and τ is the torque per unit length resulting from supplementary contributions to the internal elastic response.

Let us consider the curvilinear integral of (B.1) over each bead i, following the centerline of the skeleton joining the contact point c i-1 between bead i -1 and i and the bead center r i , as well as the bead center and the contact point c i between bead i and i + 1 (see Fig. 4). The curvilinear arclength s thus varies from 2ai to 2a(i + 1) within bead i.

Following the centerline, the integral of the internal stress contribution to (B.1) reads

(B.4)

In the limit of pointwise contacts, the normal stress produced by contact forces at contact point c i-1 located at point

where δ S stands for the surface Dirac distribution at the bead surface. Consequently the moment distribution associated with a Dirac contact forces applied at

Since the area of the cross-section S(s) normal to the centerline tends to the bead surface itself as s → 2ai or s → 2a(i + 1), one can find that n is (2ai) → f c i-1 and n is (2a(i + 1)) → f c i as a → 0 using (B.4), (B.5) and (B.2). Hence, the finite size integral of (B.1), fulfills the following limit as the bead radius tends to zero

which is consistent with the force used in [START_REF] Manghi | Propulsion with a rotating elastic nanorod[END_REF].

The second term of the moment balance equation

where volume V i-is the half-bead joining contact point c i-1 with bead center r i , whose pointing outward normal at r i

Similarly, considering the moment balance equation (B.3) from bead center r i to contact point c i leads to

where volume V i+ is the half-bead joining bead center r i to contact point c i , whose pointing outward normal at r i is -t i = e i+1,i . The surface S i+ enclosing half-bead V i+ is composed of half-sphere S i+ and disk D i+ , S i+ = S i+ ∪ D i+ . Hence the integrated contribution of the second term of the moment balance equation (B.3) is the sum of the right-hand-side of (B.8) and (B.9) which ought to be evaluated from the volume integral of the total stress over V i-∪ V i+ inside bead i. Since the stress tensor inside the beads is not known, it is possible to relate it to the applied normal force at bead surface. Using divergence theorem on any volume V , enclosed by surface S, one finds

where n is the normal pointing outward surface S, whilst introducing the tensor D S associated with the first moment contribution of the stress at surface S. If the surface S is the surface enclosing the considered bead, D S is the usual stress tensor, associated with the hydrodynamic interactions between the fluid and the bead. When considering hydrodynamic interactions, D S is usually decomposed into a symmetric tensor called stresslet and an anti-symmetric one called couplet.

Using relation (B.10) in (B.8) as well as (B.9), one finds the following four contributions 2a(i+1) [START_REF] Zhu | Drag of a flexible fiber in a 2D moving viscous fluid[END_REF] to the integration of the second term of (B.3). In the limit of bead radius tending to zero, then t i-1 → t i , so that the outward normal vector to D i-, t i-1 , tends to the opposite of the outward normal vector to D i+ . Since D i-→ D i+ , this implies in turn, that D D i-→ -D D i+ . Furthermore, since in the asymptotic limit of zero bead radius, S i-∪ S i+ → S, one finds that 2a(i+1)

where S is the bead surface here. For now, we concentrate on the contact forces contribution to (B.12). Using the contact surface force (B.5) it is then easy to compute the contact forces contribution to (B.12), 2a

2a(i)

Hence, result (B.13) is consistent with moment balance used in [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF].