A Second-Gradient Theory of Dilute Suspensions of Flexible Rods in a Newtonian Fluid

Most suspension descriptions used nowadays are based on Jeffery’s model or some phenomenological modifications of it that do not take into account size effects: the kinematics and stresses do not involve a micro-mechanical characteristic length and thus, the predicted rheological properties are necessarily independent of rod length. More sophisticated models able to enrich first-gradient kinematics as well as to activate rod-bending mechanisms are needed, in particular to explain the mild elasticity observed experimentally. In this paper we propose a second-gradient theory for dilute suspensions of flexible rods in a Newtonian fluid that is indeed able to activate rod bending and predict viscoelastic behaviour.


Introduction
Short fibers are widely used for reinforcing polymers. The resulting microstructure has been traditionally described at three different scales.
At the finest scale the microstructure can be described, when the suspension is dilute enough, by tracking a population of rods that move with the suspending fluid (assumed Newtonian) and orient depending on the velocity gradient according to Jeffery's equation [22] that relates the orientation evolution with the flow velocity field, in particular to its first gradient. In that case the motion and orientation of each fiber is assumed decoupled from the others. When the concentration increases rod-rod interactions cannot be neglected anylonger. They must be included into the force and moment balances of each rod in order to describe their motions.
When the representative population involves too many fibers, the computational efforts to track the population becomes unaffordable. Thus, the simple and well-defined physics must be sacrificed in order to derive coarser descriptions. Kinetic theory approaches [6,11,24] describe such systems at the mesoscopic scale. Their main advantage is their ability to address macroscopic systems, while keeping track of the fine physics through a number of conformational coordinates introduced for describing the microstructure and its time evolution. At this mesoscopic scale, the microstructure is defined from a distribution function that depends on physical space, time and a number of conformational coordinates, i.e. rod orientation in the case of slender bodies suspensions.
The moments of the orientation distribution function constitute a coarser description in general used in macroscopic modeling [1]. At the macroscopic scale the equations governing the time evolution of these moments usually involve closure approximations whose impact on the results is unpredictable [7,12,23]. Alternatively, macroscopic equations are carefully postulated in order to guarantee the model objectivity and its thermodynamical admissibility.
In the case of dilute suspensions of short fibers, the three scales have been extensively considered to model the associated systems without major difficulties. However, as soon as the concentration increases, difficulties appear [28]. In the semi-dilute and semi-concentrated regimes fiber-fiber interactions occur, but in general they can be accurately modeled by introducing a sort of randomizing diffusion term [15]. There is a wide literature concerning dilute and semidilute suspensions, addressing modeling [5,[17][18][19][20], flows [3,4,8,35] and rheology [27,29,30].
Thus, most of suspension descriptions used nowadays are based on Jeffery's model or some phenomenological modification of it that do not take into account size effects: the kinematics and stresses do not involve a micro-mechanical characteristic length and thus, the predicted rheological properties are necessarily independent of rod length.
Size effects disappear as soon as a constant gradient of the fluid velocity is postulated at the scale of the rod (firstgradient framework). In that case the rod kinematics consist of an affine transformation corrected in order to ensure rod inextensibility. Obviously, because of the transformation linearity, rod bending mechanisms are prevented. In [21], the deformation of an inextensible thread without bending rigidity was considered, and it was shown that in a flow with constant velocity gradient, deformations are activated as soon as a nearly straight initial configuration is assumed. In order to activate bending mechanisms, we considered in our former works first-gradient hydrodynamics acting on rods with nonstraight unloaded configurations [9,10]. In [34], the authors considered a three-bead-two-rod system equipped with an elastic spring that resists the rod differential rotation.
In this paper, we propose a second-gradient description of rod suspensions that allows activation of rod bending and thus naturally introduces elastic stresses, although it considers straight rods and predicts Jeffery rod kinematics.
In the next section, we develop first and second-gradient models for the kinematics of rods immersed in a Newtonian fluid, with the rods being modeled as dumbbells. We will prove that in both descriptions the rod kinematics are described by Jeffery's equation and that, with forces acting along the rod direction, rod bending cannot be activated. In Sect. 3, we consider an alternative approach in which hydrodynamic forces are distributed all along the rod length. In that case, the rod kinematics are again described by Jeffery's equation but, as soon as the second gradient is considered, rod bending is activated. From the microscopic description, we derive in Sect. 4 a second-gradient macroscopic model that couples the rod behaviour with the flow kinematics. In Sect. 5, we show that the proposed theory predicts a Maxwell linear viscoelastic behaviour. Section 6, we study the model predictions in the inception and cessation of simple shear flow.
Finally, in the spirit of [25], we briefly address in Sect. 7 the complex flow in a driven cavity assuming coupling between rod and flow kinematics.
We voluntarily omit a detailed numerical analysis of the proposed theory in order to limit the size of this paper. It will be reported in future publications. We thus mainly focus our analysis on the effect of flow kinematics on the rod conformation, that constitutes the main originality of our proposal.

Remark 1
In this paper, we consider the following tensor products: -if a and b are first-order tensors then the single contraction "·" reads (a · b) = a j b j (Einstein's summation convention); -ifa and b are first-order tensors then the dyadic product"⊗" reads (a ⊗ b) jk = a j b k ; -if a and b are first-order tensors then the cross product"×" reads (a × b) j = jmn a m b n (Einstein's summation convention) with jmn the components of the Levi-Civita tensor (also known as permutation tensor); -if a and b are respectively second and first-order tensors then the single contraction "·" reads (a · b) j = a jm b m (Einstein's summation convention); -if a and b are second-order tensors then the single contraction"·" reads (a · b) jk = a jm b mk (Einstein's summation convention); -if a and b are second-order tensors then the double contraction ":" reads (a : b) = a jk b k j (Einstein's summation convention); -if a and b are third-order tensors then the triple contraction "∵" reads (a ∵ b) = a jkm b mk j (Einstein's summation convention); -if a and b are fourth-order tensors then the fourth contraction "::" reads (a::b) = a jkmn b nmk j (Einstein's summation convention).

From First to Second-Gradient Descriptions of Rod Kinematics
We consider a suspending medium consisting of a Newtonian fluid of viscosity η in which there are suspended rigid and non-Brownian rods. We assume as first approximation that their presence and orientation do not affect the flow kinematics (this hypothesis will be relaxed later) that is defined by the velocity field v(x, t), with x ∈ Ω ∈ R d , d = 2, 3. The conformation of each rod of length 2L can be described from its orientation, given by the unit vector p located at the rod center of gravity G and aligned along the rod axis. The rod's mass is neglected and with it micro-scale inertial effects. The orientation evolution of an ellipsoidal particle suspended in a Newtonian fluid is described by Jeffery's equation [22], that for rods (ellipsoids with infinite aspect ratio) can be derived by considering the system illustrated in Fig. 1 consisting of a rod and two beads located at both rod ends where we assume that hydrodynamic forces apply. We assume that the forces that apply on each bead F depend on the difference of velocities between the fluid and the bead, the first one given by v 0 + ∇v · pL and the second one by v G +ṗL. Thus, force F(pL) reads: with ξ the friction coefficient, v 0 the fluid velocity at the rod center of gravity and v G the velocity of the rod center of gravity.
Obviously if F applies on the bead pL, then at the opposite bead −pL the resulting force reads By adding Eqs. (1) and (2) and enforcing the force balance (neglecting inertia effects, as the rod mass is assumed negligible), we obtain than implies v 0 = v G , that is, the rod center of gravity is moving with the fluid velocity. For alleviating the notation, we consider F = F(pL) and F(−pL) = −F. As the resulting torque must also vanish, the only possibility is that force F acts along the direction defined by p, that is F = λp, with λ ∈ R. Thus we can write Premultiplying Eq. (4) by p and taking into account that p · p = 1 and consequently p ·ṗ = 0, we get: that allows us to write from which we finally obtain the classical Jeffery equatioṅ p = ∇v · p − (∇v : (p ⊗ p)) p.
Remark 2 Because the factor ξ L appears in both sides of Eq. (6), the predicted rod kinematics does not involve size effects.
The forces applied at the rod ends pL and −pL are respectively λp and −λp, both directed along the rod direction and by construction auto-equilibrated.
With λ given by Eq. (5) the force F reads which is aligned, as expected, in the rod direction. Such a force generates the rod tension or compression depending on the sign of ∇v : (p ⊗ p). In the case of a flexible rod it can produce its extension or compression, the last being able to produce buckling.
Thus, a certain elasticity could be expected, however for standard materials the elasticity related to the axial tension is too small and compression mechanics leads to buckling and the consequent rod rupture.
Rod bending seems to be another plausible mechanism responsible for elastic effects, however forces aligned in the rod direction cannot activate rod bending. The introduction of such effects requires at least a second-gradient theory, as described in what follows.

Second-Gradient Description with Concentrated Forces at the Rod Beads
We now consider the same system but with a higher-order description of the fluid velocity field at the rod scale. Again, forces applied on each bead F depend on the difference of velocities between the fluid and the bead, the first one now including the second-order velocity gradient H according to v(pL) = v 0 + ∇v · pL + (H : (p ⊗ p)) L 2 and the second one given by v G +ṗL. Thus, the force F(pL) reads: where the third-order tensor H is defined by Obviously, if F applies on the bead pL, then at the opposite bead −pL the resulting force reads By adding Eqs. (9) and (10) and enforcing the force balance, neglecting again inertial effects, we obtain than implies v 0 − v G = −(H : (p ⊗ p))L 2 , that is, the rod center of gravity has a relative velocity with respect to the fluid at this position.

Remark 3
The fact of having obtained a non-zero relative velocity v 0 − v G = 0 does not imply the existence of a migration mechanism, as shown in "Appendix 1". See also [33].
Since the resulting torque must also vanish, the only possibility is that force F acts along p, that is F = λp, with λ ∈ R. Thus we can write which yields the same Jeffery equation that was obtained when considering the first-order velocity gradient: The forces being again aligned in the rod direction, one could infer that the second gradient does not suffice for activating bending mechanisms.
Consideration of a third-gradient description predicts rod kinematics that differ from those dictated by Jeffery's equation (see "Appendix 2"). However, the forces remain aligned along the rod direction and rod bending thus cannot be activated either.
Until now, forces were assumed applied at the rod ends (beads). However, we can distribute them all along the rod length as commonly considered when calculating particles motion by using DPD (dissipative particle dynamics) methods.
In the next section, we consider forces applied all along the rod length as was considered in [32]. We prove that as soon as a second-gradient description is retained, bending mechanism appears naturally.

Second-Gradient Description of Rods with Distributed Forces
We consider now the system illustrated in Fig. 2 consisting of a rod and the hydrodynamical forces applied all along its The resultant force F must vanish, that is that leads to the following expression for the sliding velocity (rod-fluid relative velocity at the rod center of gravity): Thus, the distributed force reads: which leads to the moment m(s) from which we can evaluate the resultant moment M with α = +L −L s 2 ds = 2 3 L 3 . This again yields Jeffery's equatioṅ i.e. the same equation as that obtained by assuming forces applied at the rod beads.
Introducing Jeffery's equation (20) into the distributed force expression (17), we obtain This force can be decomposed into two components, i.e. one, f (s), aligned with the rod and the other, f ⊥ (s), perpendicular to it: and Remark 4 We can notice from Eq. (23) that, when considering distributed forces within a first-gradient description (H = 0), the perpendicular component of the resulting distributed forces vanishes, i.e. f ⊥ (s) = 0 and bending mechanisms are once again absent. However, bending seems possible as soon as second-gradient descriptions implying H = 0 are retained.
The resultant axial force F reads expression that corresponds to the one obtained in the case of concentrated forces if we consider the following relation between the distributed and concentrated friction coefficients: The main difference with respect to the situation in which forces were assumed applied at the rod beads, occurs when considering the distributed force perpendicular to the rod. In this case in agreement with the results found when considering concentrated forces, but the distributed force implies the existence of a bending moment M(s)k (acting in the out-of-plane direction defined by the unit vector k) When the rod is rigid there are no major consequences, but in the case of flexible rods this force creates rod bending with the associated elastic effects.

Remark 5
In the case of elastic rods, the bending moment implies rod curvature. Within the small strain and displacement hypotheses, the curvature is given by where u(s) is the rod deflection with respect to its undeformed configuration, E the elastic modulus and I the moment of area of the rod cross section with respect to the out-of-plane direction. By integrating this equation twice, we can obtain the rod bent configuration that, as expected, is symmetric with respect to the rod center of gravity.
It is important to notice that moments imply an extra length L with respect to forces, and moreover each integration step for moving from curvature to displacement introduces the length L as factor.
Having proved that second-gradient descriptions activate bending effects, we propose in the next section a theoretical model able to couple flow and rod kinematics within a second-gradient framework.

Second-Gradient Flow Model of Dilute Suspensions
Involving Flexible Rods

Second-Gradient Fluid Description
Within the first-gradient formulation, the internal power for a Newtonian fluid W int reads where T = −pI + τ with τ = 2ηD. The associated balance of momentum reads: Fried and Gurtin [16] proposed a second-gradient formulation involving the vorticity gradient. They proposed a nonstandard form of the principle of virtual power with the internal power given by where ω is the vorticity vector and G is the so-called hyper-stress. The associated generalized momentum balance reads: The following constitutive equation was assumed in [25] for the fluid hyper-stress: where the parameter ι ∈ [−1, 1] controls the asymmetry of the hyper-stress and ensures a non-negative dissipation. In the previous equation L 2 f > 0 is known as the fluid gradient length.
In the case of an incompressible fluid, the introduction of both the stress and the hyper-stress constitutive equations into the generalized momentum balance leads to: where ηL 2 f represents the so-called hyper-viscosity [25]. See [25] for a discussion on the associated boundary conditions. The higher-order velocity derivatives involved in Eq. (35) require the enforcement of a larger number of boundary conditions compared to the standard first-gradient formulation.

Second-Gradient Description of a Dilute Suspension of Flexible Rods
Inspired by the use of the vorticity ω in the above secondgradient fluid flow description, we consider now the rod beads being subjected to two actions, as depicted in Fig. 3: -A hydrodynamic force F that depends on the difference of velocities between the fluid and the bead. The fluid velocity at the bead position v(pL) taking into account second-gradient effects reads The bead velocity is again given by v G +ṗL. Thus, the resulting hydrodynamic force reads -A hydrodynamic torque M that depends on the difference between the differential vorticity at the bead position and the bead rotary velocityφ. The differential vorticity is the difference between the vorticity existing at the bead position minus the one existing at the rod center of gravity. This torque could have its origin in the distributed forces applied along the rod length as illustrated in Sect. 3. Thus, the resulting torque reads where ξ R is the rotary friction coefficient.
implying a second-order relative velocity of the rod center of gravity with respect to the fluid velocity at this position. Again, for notational simplicity, we consider F = F(pL) = −F(− pL). As soon as we consider a first gradient of the vorticitybased bending mechanism ∇ω · pL, it is easy to prove (see "Appendix 3") that the rod kinematics remains unchanged relative to Jeffery's model, and that this term only affects rod bending. Thus, as proved in "Appendix 3", we obtain: -The standard expression of the forces applying on the rod beads: -Jeffery's rod kinematics: -The bending mechanism. By considering the notatioṅ ϕ =φ(pL) = −φ(− pL) and assuming a linear elastic behavior of the rod, the relation between the applied torque M and the bending angle ϕ at the rod bead is given by where E is the elastic modulus and I the moment of area of the rod cross section. Introducing the expression of the moment M (38) into Eq. (43), we obtain whose integration allows calculating ϕ and from it the moment M.

Remark 6
Tt is noteworthy to note that: -In first-gradient flows, rods orient according to Jeffery's equation without experiencing bending effects because the vorticity gradient vanishes. -In second-order flows, rods orient according to Jeffery's equation experiencing a bending induced by the second gradient. -In rigid kinematics, rod moves with the fluid without experiencing any relative movement. Because the vorticity gradient vanishes there are not bending effects, ensuring the formulation objectivity.

Rod Kinematics
The rod kinematics is then fully described by p and ϕ. Knowing p, we can locate both rod beads. Then, knowing the beads bending angle ϕ and all the forces and moments being applied at both ends, the deformed rod configuration is parabolic and symmetric with respect to the center of gravity. Thus, from (p and ϕ), we can predict the bent configuration. The evolution of both descriptors are given by Eqs. (41) and (44). For illustrating the behaviour, we consider Poiseuille's flow defined in "Appendix 1". A rod is tracked all along its pathline. In order to emphasize the representation, we make use of non-physical parameters in the orientation The integration of Eqs. (41) and (44) yields p(t) and ϕ(t), that allow us to depict the deformed rod at the position of its center of gravity x G (t). Figure 4 shows the rod evolution. In order to emphasize the configuration evolution, we superpose in Fig. 5 all rod centers of gravity.
We can notice that, as dictated by the Jeffery equation, the rod rotates counterclockwise. Moreover, the gradient of vorticity activates rod bending, mainly when that gradient is maximum θ = π 2 . Then the rod approaches the equilibrium steady-state orientation θ = π , its rotary velocity decreases and bending relaxes because the vorticity gradient vanishes progressively when approaching the flow direction.

Induced Stresses
The stress has two components, the standard one and the one related to the hyper-stress. The first one is related to forces applied at both opposite beads acting in the rod direction and includes the suspending medium contribution τ f = 2ηD: with the total stress T given by which, as expected, is symmetric. The second one is related to second-gradient fluid contribution and rod bending. We consider a partition of the total hyper-stress G consisting of the fluid G f and the rodsG r contributions, with G = G f +G r .
We consider a standard form of the second-gradient fluid contribution G f , where L 2 f is the fluid gradient length [14]. On the other hand, we compute the rod contribution to the hyper-stress G r from from which we assume the more general form

Mesoscopic Description
Making use of the orientation distribution function Ψ (x, t, p, ϕ) (see "Appendix 4" for a summary on micro-to-macro descriptions), we can proceed to the calculation of the stress within a continuous macroscopic description. When considering the stress expression (45) and the use of the quadratic closure relation A ≈ a ⊗ a, we obtain: with T = −p I + τ . On the other hand, the hyper-stress is given by: and, by using the expression of the moment, In the previous expressions, the distribution function Ψ = Ψ (x, t, p, ϕ) (its dependence on the different coordinates is omitted for the sake of clarity) gives the fraction of rods that at position x and time t have a conformation given by (p, ϕ). The domains S and C refer respectively to the domains in which coordinates p and ϕ are defined.
In order to close the formulation, we consider the secondorder tensor g (such that G r =κ g): and compute its time derivative taking into account the expressions ofṗ (41) andφ (44): If we adopt the following closure relation: then we obtain a closed form for the evolution of g: Multiplying Eq. (56) byκ, we obtain the time evolution of the contribution of rods to the hyper-stress: If we consider the coefficient affecting the vorticity gradient in Eq. (57) and take into account the relationsκ = κ and κ = E I /L, we obtain where L 2 r represents the rod gradient length. The last term in Eq. (57) involves the coefficient κ/ξ R with reciprocal time units, which, in absence of flow, controls the relaxation to the undeformed reference configuration. Thus, the inverse of this coefficient has the meaning of a relaxation time that we denote by T . Moreover, taking into account the symmetry of tensor a, we have ∇v : a = D : a. Thus, Eq. (57) can be rewritten aṡ It is important to notice that the above equation involves two closure relations, the one related to the fourth-order orientation tensor A and the one expressed in Eq. (55). The macroscopic flow model can thus be summarized as follows: -Generalized momentum balance: -Mass balance -Constitutive equation: with τ = τ f + τ r : and -Hyper-stress G = G f +G r : and witḣ

Linear Viscoelastic Behavior Predicted by the Proposed Theory
Complex fluids are tested in the linear viscoelastic regime by applying a small-amplitude oscillatory flow in order to evaluate the frequency dependence of the complex modulus. The latter has a real part, the so-called storage modulus G , and an imaginary part called loss modulus, G . Obviously, in viscous fluids and elastic solids, the storage and loss modulus vanish respectively. In general, viscoelastic materials, both behaviours coexist.
In order to perform a linear viscoelastic analysis for standard (first-gradient) fluids, it suffices to apply a smallamplitude oscillatory simple shear flow, since only firstorder derivatives of the velocity field are involved in the model. In the present second-gradient framework, however, that involves second-order derivatives, the small-amplitude oscillatory flow must contain a richer kinematics.
Thus, we enforce the following displacement field where f (y) is a function depending on the y-coordinate, that can be differentiated at least twice, and δ 0 is the oscillation amplitude. The associate velocity field is given by and its gradient reads where f (y) = d f (y) dy . The flow is applied on an initially isotropic suspension and we can ignore the flow-induced orientation (the smallamplitude deformation implies a negligible change of the orientation state). Thus, we can assume that the orientation tensor a remains isotropic at all times: where I denotes the unit tensor. The suspension constitutive equation involves a purely viscous contribution, that vanishes as soon as the flow stops and that consequently only concerns the loss modulus G (i.e. it is not involved in the elastic behavior). For this reason, in what follows, we focus on the second contribution to the stress, more precisely on the hyper-stress G, G = G f +G r , particularly on the one associated to the presence of rodsG r , because the one related to the fluid G f is purely viscous, thus vanishing as soon as the flow stops.
When focusing on the rod contribution to the hyper-stress witḣ we can neglect the first term of the right-hand side, because it only ensures objectivity of the evolution equation in the non-linear regime. Thus, the linearized problem readṡ In Eq. (76) the term (D : a) = 0 and the vorticity ω is given by which implies and consequently, We assume that the response G r,lin to the enforced oscillation has the same frequency (in view of the linearity of the model when assuming small-amplitude oscillations) but can exhibit a phase delay ϕ, that is where the complex amplitude G 0 is such that G 0 = G 0 e −iϕ . Introducing expressions (79) and (80) into Eq. (76), we obtain or with α = EL 2 r . Algebraic manipulations in Eq. (82) yield The only non-zero component of G 0 is G 0,23 , which by taking into account the expression of A 23 from Eq. (79) reads with β = − αT δ 0 f (y) 3 . When considering a parabolic viscoemtric flow, we have f (y) = cst.

Remark 8
The contribution of the hyper-stress to the tension on a plate surface with unit normal n is given by n × ∇ ·G r . For a surface defined by n T = (0, 1, 0), the only component that resists the flow v T = (u, 0, 0) isG r 32 , that results from G r 23 and G r 32 , the first of them being in the present case nonzero.
We can conclude from the above analysis that the storage modulus G = , responsible for the elastic behaviour, has a dependence on the frequency ω, exhibiting a slope of 2 at low frequencies and becoming constant (plateau) at high frequencies, i.e. the well-known Maxwellian behaviour observed in many complex fluids.

Inception of Shear Flow Followed by Relaxation
In this section, we wish to visualize from a physical point of view the subtle coupling between flow and rod elasticity, by considering two scenarios: (i) first, a prescribed shear flow acting on initially-straight flexible rods and yielding rod bending, followed by (ii) the progressive cessation of flow once the bent rods relax to their straight configurations after the flow source (pressure drop) has been abruptly set to zero.

Inception of Shear Flow
We prescribe the shear flow with f (y) at least four times differentiable since the secondgradient flow model involves the curl of the hyper-stress divergence, ∇ × ∇ · G, with the hyper-stress depending on the vorticity gradient (see Eqs. (33) and (34)). We focus on a rod located on the streamline y 0 = cst, such that f (y 0 ) = du dy | y 0 < 0.
The velocity gradient, vorticity and vorticity gradient read respectively ∇v = and We assume that at the initial time t = 0 the rods are fully contained in the x − y plane. Thus, the orientation tensor has the form where a 0,12 = a 0,21 and a 0,22 = 1 − a 0,11 . The hyper-stress evolution equation readṡ At t = 0, assuming a fully relaxed system (straight rods), G r (t = 0) = G r 0 = 0, the hyper-stress evolution only depends on the term EL 2 r a 0 · (∇ω) T , i.e.
which, taking into account Eqs. (88) and (89), giveṡ We consider two scenarios: -All rods are initially aligned in the flow direction, implying and according to Eq. (92)Ġ r 0 = 0. In this case, rods remain oriented in the flow direction (according to Jeffery's equation) and they do not experience bending because G r (t) = G r 0 = 0.
-All rods are initially aligned perpendicularly to the flow direction, implying and thus according to Eq. (92), Consequently, the hyper-stress is given at time δt by In order to interpret the physics of this result, we come back to the microscopic definition of the hyper-stress related to the rods (Eq. (52)), all of them being aligned initially in the direction given by the unit vector p 0 which implies and thus 23 . With f (y 0 ) < 0, we have G r δ,23 > 0 and then ϕ > 0 inducing rod bending as illustrated in Fig. 4. We can go one step forward, integrating Eq. (92) from t = δt to t = 2δt. Now, G r δ is given by Eq. (96). Integration of Jeffery's equation results in an updated orientation p(t = δt) = p δ given by where p δ,y ≈ 1 and p δ,x < 0, with ||p δ || = 1. Now the term a δ · (∇ω) T activates the evolution of G r 13 and G r 23 and the objectivity of the resulting evolution is ensured by the term ∇v · G r δ in Eq. (92) that takes into account that the rod is rotating counterclockwise, in order not to include the rod rotation in the bead bending angle.

Progressive Cessation of Flow Induced by Rod Relaxation
We consider now the state reached by the system at time δt with all rods initially aligned perpendicularly to the flow direction (scenario discussed in the second item of the previous section).
In that case the hyper-stress reads Neglecting G f (i.e. putting L f = 0) and for the sake of clarity neglecting also the contribution of rods to the stress (i.e. considering N p ≈ 0), the flow model (38) reduces to Thus, if we remove suddenly the external action applied for creating the flow, ie. the pressure drop ∂ p ∂ x is set to zero, the flow persists until the full relaxation of the hyper-stress, that is until rods recover their initially-straight configuration. In absence of viscosity, this relaxation occurs instantaneously, but the presence of the viscous fluid requires that flow occurs for accommodating relaxation of rod bending.
This mechanism, that can be physically fully interpreted, is due to the transient nature of the equation that governs the rod contribution to the hyper-stress. Assuming G f = 0 and G r = 0, as soon as the external action creating the flow is removed, the flow stops suddenly, because the fluid hyperstress is not concerned by a relaxation mechanism, that is, it relaxes instantaneously, like a viscous stress.

Complex Flow Simulation
In this section, we solve the proposed suspension flow model in the driven cavity flow problem. The suspension is confined in a cavity Ω = [0, L] 2 , with L = 2. The velocity is specified at the upper wall as v T (x, y = L) = (Lx − x 2 ), and it is set to zero in the remaining part of the domain boundary (left, right and bottom walls). Moreover, as proposed in [25], we enforce the extra-condition n × G · n = −ηlω × n in the whole domain boundary ∂Ω, where l andηl are respectively the so-called adherence length and the boundary viscosity. In order to take the incompressibility constraint ∇ · v = 0 into account, we make use of the stream function formulation. The stream function Ξ is related to the components of the velocity vector according to Both velocity and vorticity are written in terms of the stream function Ξ . It is easy to verify that Due to the high order of the resulting partial differential equation, a spectral Chebyshev collocation technique is used here for discretizing the motion equation.
The computed streamlines related to the second-gradient flow model solution are depicted in Fig. 6. In this simulation, we considered L f = 0, N p ≈ 0 and L r 1, which implies that the flow kinematics are very close to the ones related to a Newtonian fluid. Our main interest when considering the second-gradient description concerns the evaluation of the rod elastic loading, more than the analysis of the perturbed flow kinematics analyzed in detail in [25]. Figure 7 shows the absolute value of the rod bending angle ϕ at two different times. The initial orientation state is full alignment in the y-direction everywhere in the computational domain Ω.

Conclusions
We have developed first and second-gradient models for the kinematics of rods immersed in a Newtonian fluid, with the rods being modeled as dumbbells. It has been shown that in both descriptions the rod kinematics are described by Jeffery's equation and that, with forces acting along the rod direction, rod bending cannot be activated. We have then considered an alternative approach wherein hydrodynamic forces are distributed all along the rod length. In that case, the rod kinematics are again described by Jeffery's equation but, as soon as a second-gradient description is considered, rod bending is activated. From this microscopic description, we have derived a second-gradient macroscopic model that couples the rod behaviour with the flow kinematics. The proposed theory has been shown to predict a Maxwell linear viscoelastic behaviour. We have studied the model predictions in the inception and cessation of simple shear flow. Finally, we have illustrated the use of the proposed model in the simulation of complex flow in a driven cavity, assuming coupling between rod and flow kinematics. A deeper analysis of the flow-induced bending and the effects of rod bending on flow kinematics is work in progress. The proposed model is also being extended to concentrated suspensions involving entangled flexible rods wherein the rotary drag is created by the rotation of the rod network. This extension of the proposed theory could constitute a valuable approach for modeling complex flows of concentrated suspensions involving flexible rods as encountered in SMC processes in composite manufacturing.

Appendix: On the Relative Velocity Between the Fluid and the Rod Center of Gravity
We showed in Sect. 2.2 that, when considering a secondgradient description, the rod center of gravity has a velocity different from the local fluid velocity. Even though one could think that this relative velocity implies rod migration, we here show that it is not the case. The relative velocity derived in Sect. 2.2 is given by Let us consider a 2D Poiseuille flow described by the following kinematics: with y ∈ [−H, H ].
In this case, the only non-zero component of tensor H is: which acts on (p ⊗ p) 22 . With p T = (cos θ, sin θ), we obtain This shows that rods do not move towards regions of lower shear rates. They never leave their trajectories, y = cst, but then Eq. (120) can be rewritten as: from which we obtain the rod rotary velocity Taking into account Eq. (121) and the triple vector product property a × b × c = b(a · c) − c(a · b), we havė Now, coming back to Eq. (37) and taking Eq. (125) into account, As we can notice, the force acting on the beads has a component in the rod direction given by and another component perpendicular to it and contained in the plane where the rod orientation is defined: Since F ⊥ (pL) = −F ⊥ (−pL) and making use of symmetry considerations, equality of torques and bending angles at both rod beads are expected, that is: and respectively. Thus, we findφ =φ. This result, however, is not compatible with Thus, the only possibility is thatφ(pL) = −φ(−pL), implying F ⊥ (pL) = F ⊥ (−pL) = 0 (which yields the standard Jeffery solution for the rod rotary velocityṗ) and M(pL) = −M(−pL), which fully agrees with Eq. (131).

Kinetic Theory Description
In view of the very large number of rods involved in a typical suspension, the description based on the tracking of each individual particle, despite its conceptual simplicity, fails to address the situations usually encountered in practice. For this reason, coarser descriptions are preferred. The first plausible coarser description applies a zoom-out, in which the rod individuality is lost in favour of a distribution function.
For the sake of simplicity, we consider rigid rods in what follows. Thus, one could describe the microstructure at a certain point x and time t from the orientation distribution function Ψ (x, t, p) that gives the fraction of rods that at position x and time t are oriented in the direction p. Obviously, function Ψ verifies the normality condition: S Ψ (x, t, p) dp = 1, ∀x, ∀t, where S is the surface of the unit ball (circumference of unit radius in the 2D case and spherical surface of unit radius in the 3D case) that defines all possible rod orientations. Thus, we have substituted the fine-grain microscopic description, that requires the compuation of each individual rod orientation, with the scalar multidimensional function Ψ . However, in order to use it, one needs to derive the equation governing the evolution of the orientation distribution function Ψ .
The balance ensuring conservation of orientation probability reads: while, for inertialess rods and considering a first-gradient descriptionẋ = v(x, t), the rod rotary velocity is given by Jeffery's equation: Equation (133), known as the Fokker-Planck equation, is a well-balanced compromise between the macroscopic scale that defines the overall process, where the space and time coordinates are defined, and a finer microscopic description for the flow-induced orientation of individual rods given by the Jeffery equation (134).
The price to pay is the increase in the model dimensionality, since the orientation distribution is defined in a highdimensional domain of dimension 5 in the general 3D case, i.e. Ψ : (x, t, p) → R + where x ∈ Ω ⊂ R 3 , t ∈ I ⊂ R + , p ∈ S.
The extra-stress due to the presence of such a population of rods is determined by adding the contribution of each rod to the stress: τ = τ f + τ r = 2ηD + 2ηN p ∇v : S p ⊗ p ⊗ p ⊗ p Ψ (x, t, p) dp where A represents the fourth-order moment of the orientation distribution function and N p is the so-called particle number that depends on rod concentration and on the friction coefficient. In view of the symmetry of A, the gradient of velocities in the previous expression can be replaced by the rate of strain tensor D:

Macroscopic Description
Fokker-Planck based descriptions are rarely used in practice in view of the curse of dimensionality that the introduction of conformation coordinates implies (here, the rod orientation). Thus, standard mesh-based discretization techniques, such as finite differences, finite elements or finite volumes, fail when addressing models defined in high-dimensional spaces. For this reason, mesoscopic models are usually coarsened one step further in order to derive macroscopic models defined in standard physical domains, involving only space and time coordinates.
In this section, we illustrate the transition from the mesoscopic to the macroscopic scale. At the macroscopic scale, moments of the orientation distribution are used for describing the microstructure.
A closure relation is needed in order to express the fourthorder moment A as a function of the lower-order moments (as odd moments vanish because of the symmetry of Ψ , the only non-zero lower moment is the second-order moment a). Different closure relations have been introduced and widely used [2,13,26,31]. For example, in the quadratic closure relation (that is exact only when all rods are locally aligned in the same direction), the fourth-order moment is approximated as We thus havė a ≈ ∇v · a + a · (∇v) T − 2 (∇v : a)a, and the stress tensor is given by