Marta Perez

Adrien Scheuer

Emmanuelle Abisset-Chavanne

F Chinesta

R Keunings

A multi-scale description of orientation in simple shear flows of confined rod suspensions

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fibre suspensions can be described at different scales: (i) microscopic, the scale related to the fibre; (ii) mesoscopic, the scale related to the particle population within a representative volume of the local macroscopic conditions; and (iii) the macroscopic scale that is related to the process and the final composite part.

Suspensions involving particles can be described at the microscopic scale by tracking the motion of each individual particle involved in the system. This approach is based on two main elements: (i) the knowledge of the equation governing the particle motion in the fluid flow, and (ii) the availability of computational resources for tracking efficiently millions of particles. In dilute suspensions, the motion of ellipsoidal particles immersed in a Newtonian fluid can be accurately described by using Jeffery's equation [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] . For circumventing the difficulties related to simulations at the microscopic scale where too many particles are present, these difficulties being more computational than conceptual, coarser models were introduced.

Mesoscopic kinetic theory models result from coarsening microscopic descriptions. In kinetic theory models, the individuality of the particles is lost in favour of a statistical description that substitutes the entities by a series of conformation coordinates [START_REF] Bird | Dynamic of Polymeric Liquid[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF] . For example, when considering a suspension of rods, the mesoscopic description consists in giving the fraction of rods that at position x and time t are oriented along direction p . This information is contained in a probability distribution function -pdfwhose evolution is governed by a Fokker-Planck equation. Fokker-Planck equations being multidimensional (in time, physical space and conformation space), standard mesh-based discretization techniques fail when addressing their numerical solution. This issue is known as the "curse of dimensionality". The direct solution of Fokker-Planck equations has been made feasible in many cases with the introduction of the Proper Generalized Decomposition (PGD) approach [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] .

At the macroscopic scale, the pdf is substituted by some of its moments. Here the level of detail and the involved physics are sacrificed in favour of computational efficiency. The equations governing the time evolution of these moments usually involve closure approximations whose impact on the results must be evaluated.

In the case of dilute suspensions of short fibres in a Newtonian fluid, the three scales have been extensively considered without major difficulties to model the associated systems. Challenges appear, however, as soon as the concentration increases. In the semi-dilute and semi-concentrated regimes, fibre-fibre interactions occur, but in general they can be accurately modelled by introducing a randomizing diffusion term [START_REF] Folgar | Orientation behavior of fibers in concentrated suspensions[END_REF] . There is a wide literature on the modelling of dilute and semi-dilute suspensions, e.g. [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF][START_REF] Hand | A theory of anisotropic fluids[END_REF][START_REF] Hinch | The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles[END_REF][START_REF] Hinch | Constitutive equations in suspension mechanics. Part I[END_REF][START_REF] Hinch | Constitutive equations in suspension mechanics. Part II[END_REF][START_REF] Petrie | The rheology of fibre suspensions[END_REF] . Available models describe quite well the experimental observations.

When fibre orientation predictions are compared with experimental results in injection processes involving concentrated suspensions of short fibres, a noticeable delay in the orientation kinematics is observed. Some ad hoc modelling approaches were proposed to delay the orientation kinematics predicted with the Folgar and Tucker model [START_REF] Folgar | Orientation behavior of fibers in concentrated suspensions[END_REF] . In these models, a fluid-particle sliding mechanism is introduced or different rate equations for the eigenvalues and eigenvectors of the orientation tensor are used [START_REF] Phelps | An anisotropic rotary diffusion model for fiber orientation in short and long fiber thermoplastics[END_REF][START_REF] Wang | Prediction of fiber orientation in a rotating compressing and expanding mold[END_REF][START_REF] Wang | An objective model for slow orientation kinetics in concentrated fiber suspensions: theory and rheological evidence[END_REF] . A more physically-based approach was proposed by Ferec and co-authors in [START_REF] Ferec | Modeling fiber interactions in semiconcentrated fiber suspensions[END_REF] , where the interaction mechanisms were taken into account within a multi-scale framework.

All these studies concerned unconfined flows, despite the fact that processes of industrial interest often involve narrow gaps where fibres may have length greater than the gap and wall effects cannot be ignored. The orientation delay was indeed observed in such confined flows (e.g. [START_REF] Le | X-ray phase contrast microtomography for the analysis of the fibrous microstructure of SMC composites[END_REF] ), wherein the orientation process could probably differ from the one predicted by the standard Jeffery equation and the mesoscopic and macroscopic models derived from it.

In the present paper, we propose a multi-scale description of rod orientation in confined conditions and analyse the impact of initial conditions on rod kinematics under confinement constraints. It is well known that the orientation kinematics of fibres located near the wall deviate from Jeffery's predictions, even when the fibres never enter in contact with the wall [START_REF] Gavze | Particles in a shear flow near a solid wall: effects of non sphericity on forces and velocities[END_REF][START_REF] Jayageeth | Dynamics of short fiber suspensions in bounded shear flow[END_REF][START_REF] Mackaplowt | A numerical study of the rheological properties of suspensions of rigid, non-brownian fibres[END_REF][START_REF] Stover | The motion of rodlike particles in the pressure-driven flow between two flat plates[END_REF] . The model proposed here does not take into account these hydrodynamic confinement effects, but rather it considers physical contact between rods and gap walls through the introduction of the contact force ensuring wall impenetrability. In order to consider the rod kinematics perturbation in absence of physical contact with the wall, we should consider an extra hydrodynamic force acting on the dumbbell beads describing the rod (see below); this force would depend on the distance to the wall and the approaching velocity. The model proposed in this work aims at capturing the first-order effects of confinement. The effect of the wall proximity on the fiber kinematics constitutes a second-order effect, to be addressed in future works.

In what follows, we use the following notation together with Einstein's summation convention:

• if a and b are first-order tensors, then the single contraction "•" reads (a • b ) = a j b j ;

• if a and b are first-order tensors, then the dyadic product "⊗" reads (a ⊗ b ) jk = a j b k ;

• if a and b are respectively second and first-order tensors, then the single contraction " • " reads (a • b ) j = a jk b k ;

• if a and b are respectively third and first-order tensors, then the single contraction " • " reads (a • b ) jk = a jkm b m ;

• if a and b are second-order tensors, then the double contraction ": " reads (a : b ) = a jk b k j ;

• if a and b are respectively second and fourth-order tensors, then the double contraction ": " reads (a : b ) jk = a ml b ml jk .

Modelling confined suspensions of rods

We consider a dilute suspension of rigid, non-Brownian, highaspect-ratio fibres suspended in a Newtonian fluid of viscosity η.

The fibres are modelled as rigid rods of length 2 L . The rod orientation is given by the unit vector p located at the rod centre of gravity G and aligned along the rod axis. Inertial effects are neglected in the sequel. We assume that the presence and orientation of the rods do not affect the flow kinematics that is defined by the velocity field v T = ( ˙ γ z, 0 , 0) describing a simple shear flow. Elongation is not considered in this work because it tends to reduce and even suppress all confinement effects when its intensity with respect to shear is large enough. In the case of pure elongation, the rods align monotonically in the extension direction and can never enter in contact with the gap walls. The flow occurs in a narrow gap We consider first an individual rod whose centre of gravity is located in the mid plane z G = 0 . Thus, both rod extremities enter in contact or lose contact with the gap walls simultaneously. Without loss of generality, we consider that the unit vector p related to bead located at p L points towards the upper wall, and consequently -p points towards the opposite wall.

× [ -H, H] , with x T = (x, y ) ∈ ∈ R
As just indicated, this work only considers confined flow exhibiting constant shear rate throughout the gap thickness in order to ensure a constant shear rate along the rod length. This flow is conceptually of interest because it allows one to identify the differences between the standard (unconfined) Jeffery model and the one that results when confinement effects take place. It remains, however, quite far from the applicative processing conditions that in general involve more complex flows, in many cases of Poiseuille type. In processing conditions, almost-parabolic velocity profiles throughout the gap thickness are usually encountered, implying that the velocity gradient is no longer constant along the length of the rod when its extremities approach both walls. In these circumstances, one should consider second-order kinematics in the derivation of the rod rotary velocity, as done in [START_REF] Abisset-Chavanne | A second-gradient theory of dilute suspensions of flexible rods in a newtonian fluid[END_REF][START_REF] Borzacchiello | Orientation kinematics of short fibres in a second-order viscoelastic fluid[END_REF] for unconfined flows. In that applicative perspective, one should also take into account non-Newtonian rheology of the suspending fluid, which adds a major difficulty as discussed in [START_REF] Borzacchiello | Orientation kinematics of short fibres in a second-order viscoelastic fluid[END_REF] . All these issues will be presented and discussed in ongoing publications. In the present work, we concentrate on simple shear flows, Newtonian suspending fluids, and we assume incompressible and isothermal conditions.

Microscopic description of confined kinematics of an individual rod

It is well known that the kinematics of a rod with infinite aspect ratio immersed in an unconfined simple shear flow with velocity gradient ∇v is given by Jeffery's equation [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] ˙

p = ∇v • p -(∇v : (p ⊗ p )) p . (1)
Jeffery's equation predicts full alignment in the flow direction. For finite aspect ratio ellipsoidal particles, periodic trajectories known as Jeffery's orbits are predicted instead. Brownian effects also avoid full alignment and are generally introduced at the rod population level.

When the domain thickness is smaller than the rod length ( H < L ), some Jeffery's trajectories are forbidden, i.e. those trajectories involving L p ( t ) • n > H , where n = (0 , 0 , 1) T is the unit vector defining the thickness direction. In that case, the rod kinematics are defined by the standard Jeffery model (1) while p • n < H / L and are perturbed as soon as the rod reaches the upper and lower walls. In order to determine the perturbed rod kinematics, we represent the rod as a dumbbell [START_REF] Abisset-Chavanne | On the multiscale description of dilute suspensions of non-brownian rigid clusters composed of rods[END_REF] with hydrodynamic and contact forces acting on the dumbbell beads as illustrated in Fig. 1 . We assume that hydrodynamic forces applied on each bead F H depend

F H (p L ) = ξ (v 0 + ∇v • p L -v G -˙ p L ) , ( 2
)
where ξ is the friction coefficient, v 0 the fluid velocity at the rod centre of gravity, and v G the velocity of the centre of gravity. The contact force is assumed to act in the direction perpendicular to the wall (we thus ignore possible friction at the wall):

F C (p L ) = μn , (3)
with

F C (p L ) = -F C (-p L )
. The contact force intensity μ is unknown at this stage and will be determined such as to prevent the beads from penetrating the walls.

Thus, the total force acting on the bead located at p L is

F (p L ) = F H (p L ) + F C (p L ) .
As inertia is neglected, the resultant force acting on the dumbbell vanishes. Since

F C (p L ) = -F C (-p L ) , the force bal- ance yields F H (p L ) = -F H (-p L ) and thus v G = v 0 , (4)
i.e. the velocity of the rod centre of gravity coincides with the velocity of the fluid at that position. This ensures that both beads are simultaneously in contact with the walls.

The resulting torque must also vanish, that is

F = F H + F C = λp , with λ ∈ R . Since v 0 = v G , we have ξ (∇v • p L -˙ p L ) + μn = λp . ( 5
)
Premultiplying Eq. ( 5) by p and taking into account that p • p = 1 and ˙ p • p = 0 , we obtain

λ = ξ L (∇v : (p ⊗ p )) + μp z , (6)
with p z = p • n . Injecting this expression for λ into Eq. ( 5) yields

ξ L (∇v • p -˙ p ) + μn = ξ L (∇v : (p ⊗ p )) p + μp z p , (7)
or

(∇v • p -˙ p ) + μ ξ L n = (∇v : (p ⊗ p )) p + μ ξ L p z p . ( 8
)
The rotary velocity ˙ p is thus given by ˙

p = ∇v • p -(∇v : (p ⊗ p )) p + μ ξ L (n -p z p ) = ˙ p J + ˙ p C . ( 9
)
Here, ˙ p J denotes the Jeffery rotary velocity component given by Eq. (1) , while the confined component is defined as ˙

p C = μ ξ L (n - p z p ) .
We now need to determine the contact force intensity μ. In order to obtain its value we must consider that the contact force appears in order to avoid that the bead leaves the flow domain, that is ˙

p • n ≤ 0 , with μ = 0 if p z L = H and ˙ p J • n > 0 .
The contact force μn must ensure that the resulting velocity is tangent to the wall. Because dynamical effects are neglected, the bead cannot rebound. It thus suffices to enforce the condition ˙ p • n = 0 . Multiplying Eq. ( 9) by n yields or

0 = ˙ p J • n + μ ξ L (1 -p 2 z ) = ˙ p J z + μ ξ L (1 -p 2 z ) , (10)
μ = - ξ L 1 -p 2 z ˙ p J z . ( 11
)
Eq. ( 9) thus reduces to

˙ p = ˙ p J - 1 1 -p 2 z ˙ p J z ( n -p z p ) = ˙ p J + ˙ p C , ( 12
)
with ˙ p J given by Eq. ( 1) ,

p z = p • n and ˙ p J z = ˙ p J • n .
The first term in Eq. ( 12) corresponds to the standard, unconfined Jeffery kinematics, while the second term avoids that rod beads leave the flow domain. While preparing this paper, it came to our attention that Eq. ( 12) had been derived independently in [START_REF] Ozolins | Simple models for wall effect in fiber suspension flows[END_REF] using similar arguments but a different procedure.

Eq. ( 9) can be rewritten as

˙ p -• p = D • p -(D : (p ⊗ p )) p + μ ξ L (n -p z p ) , ( 13
)
from which we can conclude that the rate of strain due to the flow kinematics and the reaction forces induce an effective rotation ˙ p -• p , thus ensuring the model objectivity. So, to sum up, the kinematics of a rigid rod of length 2 L in a simple shear flow occurring in a narrow gap of thickness 2 H are given by the following confined Jeffery model :

⎧ ⎨ ⎩ ˙ p = ˙ p J if p z L < H ˙ p = ˙ p J if p z L = H & ˙ p J • n ≤ 0 ˙ p = ˙ p J + ˙ p C if p z L = H & ˙ p J • n > 0 , ( 14
)
where the rods are in contact with the walls if p z = H/L . It is important to notice that expression [START_REF] Doi | The Theory of Polymer Dynamics[END_REF] only applies if the trial Jeffery velocity is such that ˙ p J • n > 0 . When p z L < H , ˙ p reduces to the standard Jeffery contribution ˙ p J . When p z L = H and ˙ p J • n ≤ 0 , the rotary velocity also reduces to the Jeffery contribution. The contact being unilateral, a rod can detach from the wall as soon as its velocity induces the detachment. These conditions, known as Kuhn and Tucker or Signorini conditions, are similar to those encountered when describing unilateral contact or elastoplasticity. When projecting Eq. ( 12) in the thickness direction, we have ˙ p J • n = 0 , thus ensuring that the rod beads do not leave the flow domain.

In the numerical tests discussed below, we will occasionally assume fully-confined conditions , meaning that fibres are assumed to be and always remain in contact with the walls, while their orientation evolution is governed by Eq. ( 12) . In such case, the contact is bilateral and fibres are prevented to detach from the wall even when ˙ p J • n ≤ 0 . A reaction force must appear to maintain contact . Confined Jeffery solution a disc ( t ) versus fully-confined macroscopic orientation tensor a C ( t ) in the case of rods aligned in the same direction. It can be noticed that the planar components of the orientation tensor evolve smoothly whereas the out-of-plane component shows a localized evolution when, after a period of time during which the fibre beads slide at the wall, they detach from it. and avoid detachment of the rod bead from the wall. In the sequel, this particular model will be referred to as the fully-confined Jeffery model . Although not entirely physical, it will be useful for validation purposes.

Mesoscopic description of a population of rods

Having described the kinematics of individual rods in confined flows, we now turn to a population of non-interacting rods. There are two natural approaches for doing so, i.e. discrete and continuous.

Discrete description

The discrete approach consists in computing the orientation of each individual rod belonging to a large discrete ensemble of N rods. Thus, the population is described from the individuals composing it, whose conformation is given by vectors p i , i = 1 , ..., N , each governed by Eq. ( 14) . The main drawback of this approach lies in the necessity of tracking the evolution of each rod by solving the corresponding equation, and even if conceptually there is no major difficulty, the computing cost could be excessive in most practical applications.

Continuous description

The continuous approach uses the pdf ψ( x , t , p ) that gives the fraction of rods that are oriented along direction p at position x and time t . This description avoids the complexity related to the immense number of fibres involved in suspensions of practical interest.

The pdf satisfies the normalization condition:

S ψ (x , t, p ) dp = 1 , ∀ x , ∀ t ( 15
)
where S is the rod conformation space, i.e. the surface of the unit sphere in the unconfined case. Conservation of probability leads to the so-called Fokker-Planck equation

∂ψ ∂t + ∇ x • ( ˙ x ψ ) + ∇ p • ( ˙ p ψ ) = 0 , ( 16
)
where the rod rotary velocity ˙ p is given by the confined Jeffery's equation ( 14) and ˙

x = v (x , t ) .
In the confined case, the permitted orientation domain is obtained by removing from the surface of the unit sphere S both polar regions located beyond the parallels z = ±H. Moreover, if the pdf is defined in the resulting 2 D surface, boundary layers are expected on parallels z = ±H, where rods orient while keeping contact with the upper and lower walls, which implies a mesh dependence of the discrete solution because rods concentrate on the parallels that have a null measure on the unit sphere. Appendix A discusses a possible alternative route for accomplishing mesoscopic modelling and simulation, however that route is not considered in the present work.

Macroscopic description

A macroscopic model describes the suspension microstructure via suitable moments of the pdf, defined in standard physical domains involving space and time. These moments can be computed either via a discrete or continuous approach.

Discrete approach

Consider for example the second-order moment of the pdf, also known as the second-order orientation tensor. In the discrete approach, we integrate Eq. ( 14) for the N rods of the population and compute at each instant the associated orientation tensor according to the ensemble average

a disc (t) = 1 N N i =1 p i (t) ⊗ p i (t) . (17)
Here, the superscript disc refers to the discrete approach.

As mentioned before, the main disadvantage of this approach is the computational cost due to the extremely large number of particles to be considered.

Continuous approach: the two limiting cases

In the continuous approach, the orientation distribution function is substituted by its moments for describing the microstructure [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF] at the macroscale, and an evolution equation for these moments is derived from the Fokker-Planck equation. Usually, macroscopic descriptions of rod suspensions are based on the use of the first two non-vanishing moments, i.e. the second and fourth-order moments, a and A , respectively defined by a = S p ⊗ p ψ (x , t, p ) dp , [START_REF] Hinch | The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles[END_REF] and Odd moments vanish in view of the symmetry of the pdf:

A = S p ⊗ p ⊗ p ⊗ p ψ (x , t, p ) dp . (19)
ψ (x , t, p ) = ψ (x , t, -p ) .
We consider two extreme situations, namely the unconfined case which yields the standard macroscopic description based on Jeffery's equation, and the fully-confined case wherein all rods are in contact with the gap walls. The latter is certainly relevant in conditions of intense confinement, namely for H L . In this case, we can further assume that rod orientations are distributed in the allowed region of the unit sphere but under a lubrication kinemat-ical constraint. These three modelling frameworks are described below.

• The unconfined case: For unconfined rods, the orientation is defined on the surface of the unit sphere S, and the associated Fokker-Planck equation is readily exploited to derive the classical evolution equation for the second-order orientation tensor associated to the standard Jeffery model [START_REF] Binetruy | Flows in Polymers, Reinforced Polymers and Composites: A Multiscale Approach[END_REF] :

˙

a J = ∇v • a J + a J • (∇v ) T -2 A J : ∇v . ( 20
)
Here, the superscript J indicates that we consider the orientation tensors associated with the standard Jeffery model for unconfined systems. • The fully-confined case: In the other limiting case, we assume that all rods in the suspension are in contact with the walls and that the contact is bilateral, i.e. rods can orient but always remain in contact with the walls. In this situation, we can derive an evolution equation for what we call the fully-confined orientation tensor a C defined as

a C = C p ⊗ p ψ (p ) dp , ( 21
) with C = { p , p z = ±H/L } . Indeed, we have ˙ a C = C p ⊗ p ˙ ψ (p ) dp , ( 22
)
where

˙ ψ = -∇ p • ( ˙ p ψ (p )) . ( 23
)
Integration by parts with respect to coordinates p yields

˙ a C = C ( ˙ p ⊗ p + p ⊗ ˙ p ) ψ (p ) dp . ( 24
)
In the confined case with p z = ±H/L, from which Eq. ( 21) gives

p = q p z , ( 25
)
a C = C q ⊗ q p z q p z q T p 2 z ψ (p ) dp . ( 26
)
Considering now Eq. ( 24) , with

˙ p ⊗ p = ˙ q 0 ⊗ q p z , (27)
and

p ⊗ ˙ p = q p z ⊗ ˙ q 0 , (28)
we obtain

˙ a C = C ( ˙ q ⊗ q + q ⊗ ˙ q ) p z ˙ q p z ˙ q T 0 ψ (p ) dp . ( 29
)
Now, we define the in-plane second-order orientation tensor b according to b = C q ⊗ q ψ (p ) dp . [START_REF] Stover | The motion of rodlike particles in the pressure-driven flow between two flat plates[END_REF] This new tensor does not have a unit trace since q 2 = 1 -H 2 L 2 . We also define a vector c , c = C q p z ψ (p ) dp , [START_REF] Wang | Prediction of fiber orientation in a rotating compressing and expanding mold[END_REF] which does not vanish because the pdf is symmetric, i.e. ψ (p ) = ψ (-p ) and q and p z have opposite sign on each parallel z = ±H defining C. With these definitions and Eq. ( 26) , we obtain the orientation tensor a C ,

a C = b c c T H 2 L 2 , ( 32
)
and its time derivative reads

˙ a C = ˙ b ˙ c ˙ c T 0 , ( 33
)
with ˙ b and ˙ c given by Eq. ( 29) . In order to obtain ˙ b and ˙ c , we need to derive the expression of ˙ q . For that purpose, we decompose the velocity gradient according to

∇v = G g j T G , ( 34
) such that ∇v • p = G g j T G q p z = G • q + g p z j T • q + G p z . ( 35
)
As detailed in Appendix B , we obtain ˙ q = G • q -δ 1 (G : (q ⊗ q )) q + δ 2 g δ 3 (q T • g ) q ,

(36)

with δ 1 = 1 1 -p 2 z , δ 2 = p z and δ 3 = p z 1 -p 2 z = δ 1 δ 2 .
Thus, the time derivative of tensor b is given by

˙ b = G • b + b • G T -2 δ 1 G : B + (g ⊗ c + c ⊗ g ) -2 δ 3 B • g , ( 37
)
with B and B being respectively the third and fourth-order in-plane orientation tensors related to the in-plane orientation q . Both will be expressed from c and b through adequate closure relations. The simplest closure consists in assuming

B = 1 p z b ⊗ c B = b ⊗ b . ( 38
)
It is exact when the rods are fully aligned (i.e. when the pdf reduces to a Dirac delta distribution). In the general case, its validity must be checked carefully. The development of more accurate closure approximations following the rationale considered in [START_REF] Chung | Improved model of orthotropic closure approximation for flow induced fiber orientation[END_REF][START_REF] Cintra | Orthotropic closure approximations for flowinduced fiber orientation[END_REF][START_REF] Kroger | Consistent closure schemes for statistical models of anisotropic fluids[END_REF][START_REF] Pruliere | Recirculating flows involving short fiber suspensions: numerical difficulties and efficient advanced micro-macro solvers[END_REF] constitutes a work in progress that will be reported in future publications. Finally, we obtain the time evolution of the first-order moment

˙ c = C ˙ q p z ψ (p ) dp , (39)
by considering the expression of ˙ q given by Eq. ( 36) ,

˙ c = G • c -δ 1 δ 2 G : B + δ 2 2 g -δ 3 δ 2 b • g . ( 40
)
• The lubrication simplified model: Inspired by lubrication theory that successfully reduces the 3D flow equations in the case of thin gaps by neglecting out-of-plane velocities, in the case of intense confinement, i.e. H / L < 0.3, we could assume ˙ p z ≈ 0 . Moreover, the in-plane components of ˙ p J and of its confined counterpart ˙ p given by ( 12) are very close in such situation. Thus, one could ignore the unconfined motion of rods in the thin gap and consider that all rods are each fully confined by an imaginary wall located at the initial bead z -coordinates. This is equivalent to computing the fully-confined solution a C at each z , a C ( t ; z ), z ∈ [0, H ], and then considering the average ˜

a C (t; H) = H 0 ω(z) a C (t; z) dz, ( 41
)
where ω( z ) is the fraction of rods that initially have a bead located at z (the other bead being at -z).

Thus, ˜ a C could be retained as a simplified orientation descriptor when confinement becomes dominant.

In what follows, ˜ a C is calculated by integrating a C at different coordinates z i uniformly distributed in the gap [0, H ] and then averaging these solutions. For an isotropic initial distribution in the gap, we obtain the through-the-thickness average of a C . Another possibility, not exploited in the present paper, consists in deriving the equation governing the time evolution of tensor ˜ a C by introducing Eqs. (37) and (40) into the time derivative of Eq. ( 41) and then performing integration. The derivation is detailed in Appendix C , where we show that, even though a closed evolution equation does not exist, suitable approximate expressions can be obtained. Comparison of confined a disc and unconfined a J , disc solutions computed by using a population of rods large enough to ensure their use as reference solutions, and the orientation tensor a J ( t ) obtained by integrating Eq. ( 20) using the quadratic and hybrid closure relations, a quad and a hybr respectively, for different initial confined orientation distributions.

Model predictions in simple shear flow

In this section, we discuss the predictions of the proposed models for the case of a simple shear flow with velocity v T = ( ˙ γ z, 0 , 0 )

and ˙ γ = 1 . Since we consider a unit shear rate, the time coordinate in most of the graphical representations that follow can be viewed as a shear strain coordinate. In all simulations carried out and discussed below, the coupling between flow and orientation is voluntarily neglected. The reasons are threefold. First, we wish to focus on the orientation process for a given unperturbed velocity field in order to analyze the confinement effects without having other disturbances than the ones related to the orientation mechanisms. Second, before addressing the semi-dilute or semi-concentrated flow regimes, we believe that the dilute case must be understood beforehand. Third, the constitutive equation relating the extra-stress to the orientation description should be revisited and probably modified in view of confinement effects. This analysis is currently underway within our group.

Evaluating the trajectory of a single rod

First, we consider the evolution of a single rod, or equivalently, of a population of rods all of them aligned in the same direction, using both the proposed extension of Jeffery's equation to confined systems and the equations governing the evolution of the moments of the distribution function under confinement conditions.

Solution of the confined Jeffery equation

Here, we track the orientation of a rod initially unconfined, that is p z L < H , but whose Jeffery trajectory implies at a certain instant that p J z (t) L > H. It is important to note that in our model all initial orientations associated with Jeffery's trajectories never reaching the domain wall (or reaching it at the highest point where the trajectory becomes tangent to the wall) will follow a Jeffery trajectory without any perturbation. As indicated in the introduction, hydrodynamic effects appear when the rod approaches the gap walls. These second-order effects are ignored in the simulations that follow.

Fig. 2 shows the orientation trajectory followed by the rod. At the beginning, because the rod is unconfined, it follows the trajectory dictated by Jeffery's equation [START_REF] Abisset-Chavanne | On the multiscale description of dilute suspensions of non-brownian rigid clusters composed of rods[END_REF] . The actual trajectory p ( t ) is depicted in red whereas the unconfined Jeffery motion p J ( t ) is represented by the broken blue curve. Until reaching the walls, as expected, both trajectories superpose and consequently cannot be distinguished. As soon as one of the rod beads reaches the upper wall p z L = H (the other bead reaches simultaneously the lower wall, but by using symmetry arguments we only refer to the one touching the upper wall), with ˙ p J • n > 0 for both trajectories, the actual trajectory and the one associated to the unconfined Jeffery motion bifurcate from each other. The Jeffery trajectory continues its unconfined motion with p J z L > H whereas the confined rod continues to slide on the upper wall p z L = H until it detaches from it.

The orientations satisfying the contact condition p z L = H define the two thick black parallels depicted in Fig. 2 . Thus, we notice that during a time interval the rod orients toward the flow direction while remaining in contact with the walls. In fact, during this period of time, the confined rod is crossing an infinity of forbidden Jeffery trajectories with ˙ p J • n > 0 . Finally, at a certain instant the rod reaches an unconfined Jeffery orbit that is tangent to the wall. This means that even if p z L = H, we have ˙ p J • n = 0 and then the contact force vanishes. The rod having thus reached a permitted, unconfined Jeffery trajectory, it follows it for the remainder of the simulation. Although the final orientation, i.e. full alignment of the rod with the flow, is the same as the one reached by the purely unconfined Jeffery motion, the confined trajectory has a slightly higher time of flight.

In order to evaluate the effect of confinement on the orientation time, we consider different confinement ratios H / L and calculate the time elapsed between the instant at which the confined and unconfined trajectories (having the same starting point) diverge (after reaching the wall) and the instant at which both trajectories reach an orientation degree quantified by p x = 0 . 9 . The different unconfined versus confined trajectories followed by the fibres are depicted in Fig. 3 , and the delay time normalized by the unconfined flight time (given by Jeffery's solution) for the different confinement ratios are reported in Table 1 .

Fully-confined orientation tensor

Here, the initial confined orientation is given by p T (t = 0) = (-p x , p y , H/L ) , with p y ≈ 0, H/L = 0 . 8 and such that p (t = 0) = 1 . The confined Jeffery equation ( 14) was integrated with the initial condition p (t = 0) . The orientation tensor at each time step a disc ( t ) was calculated from a disc (t) = p (t) ⊗ p (t) (having all rods oriented in the same direction is equivalent to consider N = 1 in Eq. ( 17) ), and compared with the fully-confined orientation tensor a C ( t ) obtained by integrating Eq. (33) from the initial condition

a C (t = 0) = p (t = 0) ⊗ p (t = 0) that allowed to define b (t = 0) and c (t = 0) .
Fig. 4 compares both solutions. It can be noticed that they are in perfect agreement at the beginning of the orientation process when confinement is intense, while they differ slightly from each other at the end when fibres detach. As just mentioned, the confined macroscopic model involving a C prevents detachment.

It is important to notice that in this case the quadratic closures (38) involved in the evolution equation (33) for a C are exact. Thus, any difference between both solutions is probably due to the fact that the description given by a C assumes fully-confined rods (i.e. they never detach from the wall), while Eq. ( 14) allows rods to detach from the wall.

In order to confirm this hypothesis on the origin of the noticed deviations, we solved the same problem but, when integrating the confined Jeffery equation, fibres were not allowed to detach from the wall and remained in contact with the wall for the remainder of the simulation. Fig. 5 shows, that when full-confinement conditions apply, both solutions are in perfect agreement. It is important also to notice that when the orientation distribution is fully concentrated, i.e. given by a Dirac delta distribution, the third and fourth-order closures of B and B respectively, involved in Eqs. (37) and (40) , are exact and then the discrete and continuous solution procedures remain in perfect agreement.

Evolution of a population of rods

We now consider a population of rods having different initial orientations and thus following different trajectories.

Discrete calculation of the orientation tensor

We consider a population of N = 40 0 0 rods having an almost isotropic initial orientation on the allowed part of the unit ball surface, the one limited by the parallels z = ±H, as shown in Fig. 6 . The trajectory of each p i ( t ) is obtained by integrating the confined Jeffery model [START_REF] Ferec | Modeling fiber interactions in semiconcentrated fiber suspensions[END_REF] . The results will be compared with those of the unconfined case, p J i (t) , governed by the standard Jeffery equation [START_REF] Abisset-Chavanne | On the multiscale description of dilute suspensions of non-brownian rigid clusters composed of rods[END_REF] .

At each instant, the second-order orientation tensor is computed for both the confined and the unconfined systems, a disc ( t ) and a J , disc ( t ) respectively, according to:

a disc (t) = 1 N N i =1 p i (t) ⊗ p i (t) , (42) and a J,disc (t) = 1 N N i =1 p J i (t) ⊗ p J i (t) . ( 43
)
In Fig. When considering a population of rods rather than a single rod, the delaying effect of confinement is less noticeable due to the averaging of results involving a number of unconfined trajectories associated with rods that never reach the walls. Thus, we can conclude that the impact of confinement is quite moderate in what concerns the orientation moments.

Impact of the closure relation

We now assume a population of rods with a fully-confined initial isotropic distribution, isotropic in what concerns the in-plane confined tensor b (t = 0) as depicted in Fig. 8 .

The solution of the fully-confined macroscopic model for a C ( t ) is compared with the integration of the confined Jeffery model [START_REF] Ferec | Modeling fiber interactions in semiconcentrated fiber suspensions[END_REF] . Fig. 9 compares both solutions and differences are noticed from the very beginning of the orientation process. These differences could be attributed to the closure relations (38) involved in the formulation of a C . It is easy to verify that the closure relations previously introduced for expressing B and B as a function of b and c are only exact in the case of full alignment (when the orientation pdf reduces to a Dirac distribution).

In order to quantify more precisely the impact of the closures on the computed solution a C , Fig. 10 compares a C with the secondorder orientation tensor obtained from the integration of the fullyconfined Jeffery equation ( 12) for all fibres in the population, integration that prevents their detachment from the wall. From these results, we conclude that errors introduced by the closure relations (38) accelerate noticeably the dynamics of alignment in confined suspensions, but cannot explain entirely the differences that Fig. 9 reveals.

Unconfined versus confined Jeffery models

The results just discussed indicate that confinement only has a slight influence on the kinematics of the orientation process, as Fig. 7 reveals. It is thus of interest to have a further look at the unconfined Jeffery equation (1) and the associated evolution equation [START_REF] Hinch | Constitutive equations in suspension mechanics. Part II[END_REF] for the orientation tensor a J when considering a confined initial orientation distribution. Here, we use the quadratic and hybrid closure relations for expressing the fourth-order orientation tensor A J in terms of a J [START_REF] Petrie | The rheology of fibre suspensions[END_REF] . Other closures exist, such as the natural and the orthotropic ones [START_REF] Chung | Improved model of orthotropic closure approximation for flow induced fiber orientation[END_REF][START_REF] Cintra | Orthotropic closure approximations for flowinduced fiber orientation[END_REF][START_REF] Dupret | Modelling the flow of fibre suspensions in narrow gaps[END_REF] . The former does not, however, have an explicit expression in the 3D case and the latter predicts spurious oscillations in absence of diffusion. For these reasons, we consider the quadratic are hybrid closures that are expected to perform reasonably well in unconfined conditions.

In the numerical experiments that follow, we specify initially two Gaussian orientation states on the allowed part of the unit ball surface, both depicted in Fig. 11 , limited by the parallels z = ±H, with H/L = 0 . 2 and 0.8. Thus, for integrating Eq. ( 20) , we specify

a J (t = 0) = 1 N N i =1 p i (0) ⊗ p i (0) , (44)
where the initially-confined discrete orientations p i (0) are those shown in Fig. 11 .

Fig. 11 compares the confined and unconfined solutions computed by using a population of rods large enough to ensure their use as reference solutions, and the orientation tensor a J ( t ) obtained by integrating Eq. ( 20) using the quadratic and hybrid closure relations, for the two initial confined orientation distributions depicted in that figure. The main conclusions that can be drawn from these results are:

• The out-of-plane component ( • zz ) for the unconfined Jeffery model is only slightly different from the one related to the confined model; • For confined initial conditions, confined kinematics delay slightly the evolution of the second-order orientation tensor obtained from a population of rods large enough to consider the computed solution as almost exact statistically (i.e. solving the corresponding Fokker-Planck equation would give the same result). The delay increases with the degree of confinement, as measured by the reciprocal of the ratio H / L ; • Closure relations are responsible for artificially accelerating the evolution of the second-order moment resulting from the integration of its unconfined evolution equation with confined initial conditions. Closure approximations become less accurate (i.e. the orientation process accelerates) as the degree of confinement increases (i.e. as H / L decreases); • The quadratic and hybrid closures produce similar results.

Two main reasons could be advanced for explaining the noticed deviations: (i) standard closure relations are inappropriate in presence of confinement; and (ii) the second-order moment alone is not an accurate descriptor of a constrained probability distribution.

A natural route is the development of empirical closure relations, inspired from the works of [START_REF] Chung | Improved model of orthotropic closure approximation for flow induced fiber orientation[END_REF][START_REF] Cintra | Orthotropic closure approximations for flowinduced fiber orientation[END_REF][START_REF] Dupret | Modelling the flow of fibre suspensions in narrow gaps[END_REF]] that were successfully implemented in [START_REF] Kroger | Consistent closure schemes for statistical models of anisotropic fluids[END_REF][START_REF] Pruliere | Recirculating flows involving short fiber suspensions: numerical difficulties and efficient advanced micro-macro solvers[END_REF] . This could be a valuable route indeed, and a fitted closure should of course work well as long as the operating conditions remain similar to the ones that served to construct it. Although useful from an applicative point of view, this approach has two main limitations: (i) it could hide the real physical reasons for the noticed deviations; and (ii) a closure fitted from a particular confined flow could fail as soon as the flow and degree of confinement differ significantly from the ones that were used to obtain the particular fit.

In order to analyze the ability of a second-order tensor to represent a confined orientation distribution, we decided to compute the solution of a J ( t ) from Eq. ( 20) , but instead of using a closure relation for the fourth-order orientation tensor A J , we obtained it from the solution of the equation governing its time evolution that in its turn involves the sixth-order orientation tensor A J (t) . In order to avoid once again the use of a closure relation, we computed A J (t) from its evolution equation that now involves the eightorder orientation tensor. We finally substituted the latter with the closure approximation A J (t) ⊗ a J (t) . Despite these effort s to consider a formulation up to order 6, involving a J , A J and A J with their 3 2 + 3 4 + 3 6 components and their corresponding evolution equations (in fact much less, due to normalization conditions and symmetry properties), the computed solution was found to be very poor.

Following the rationale described in [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF] , we attempted to describe the initial confined distribution depicted in Fig. 6 from its associated moments a , A and A . Use of the first three nonvanishing moments was found insufficient for approximating the actual distribution. Many additional higher-order moments would be needed, which is not practical and would compromise the efficiency of the solution procedure. As discussed in [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF] , the orientation tensors can be viewed as the Fourier series expansion coefficients of the orientation distribution function. Limitations of Fourier series for approximating rectangular functions like the ones associated to a uniform distribution in a region of the conformation space (surface of unit sphere) are well known. In that situation, it is preferable to consider the direct solution of the Fokker-Planck equation, which constitutes an appealing route for further developments [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] .

Lubrication approach

In order to validate the lubrication approach of the fully confined macroscopic model ( Section 2.3.2 ), we compare in Fig. 12 for two different degrees of confinement, H/L = 0 . 2 and H/L = 0 . 3 , the solutions obtained from the confined Jeffery model ( 14) that is considered as the reference solution, associated to an initially uniform orientation distribution in the allowed region of the unit sphere ( Fig. 6 ), and the lubrication counterpart of the fullyconfined macroscopic model. It can be noticed that the orientation kinematics is described quite well by the proposed macroscopic model, which thus seems a valuable tool for calculating the orientation evolution in highly-confined systems.

In order to understand the lubrication mechanisms and more concretely the effects of averaging between different fully-confined orientation evolutions, we depict in Fig. 13 the fully-confined solutions obtained for different confinement ratios H / L .

Conclusions

In this paper, we have extended the standard Jeffery model for rod kinematics in a Newtonian fluid in order to take account of confinement effects in a simple shear flow occurring in a narrow gap. The proposed confined Jeffery model ( 14) is meant to describe the kinematics of individual rods in confined flows. It is in principle easily exploited for a discrete description of a population of non-interacting, confined rods ( Section 2.2.1) . The development of a continuous mesoscopic description is more delicate. Although a sophisticated, confined Fokker-Planck model, based on the use of two probability distribution functions, is derived in Appendix A , we did not pursue this route in view of the anticipated numerical challenges. Finally, we have developed a continuous model for the macroscopic scale in the limiting case where rods can orient but always remain in contact with the walls. Using suitable closure approximations, we have thus obtained the evolution equation (33) for the fully-confined second-order orientation tensor a C , as detailed in Section 2.3.2 and Appendix B . In the case of intense confinement, i.e. H / L < 0.3, we can ignore the unconfined motion of rods in the thin gap and consider, in the spirit of lubrication theory, that each rod is fully confined by an imaginary wall located at its extremities. This led us to retain the gap-averaged, fullyconfined orientation tensor ˜ a C as suitable descriptor when confinement is dominant. It can be computed either via the average (41) involving values of a C through the gap, or else as solution of an evolution equation derived in Appendix C .

Numerical experiments have been conducted for a simple shear flow, with rods having initial orientations on the allowed part of the unit ball surface, i.e. the one limited by the walls. We found that the orientation kinematics predicted with the confined and unconfined Jeffery models applied to a population of rods are quite similar, with only a slight delay in confined systems.

In view of this result, we performed the same simulations with the unconfined Jeffery equation ( 1) and the associated macroscopic evolution equation [START_REF] Hinch | Constitutive equations in suspension mechanics. Part II[END_REF] for the second-order orientation tensor a J .

To our surprise, radically different results were obtained: the evolution of orientation as predicted by the macroscopic model for the orientation tensor a J ( t ), i.e. the second moment of the orientation pdf, is much faster than that obtained by computing with the Jeffery model the orientation evolution a J , disc ( t ) of a discrete population of rods. The origin of this difference is due to the impossibility of describing the confined orientation pdf and its time evolution using only the second-order moment of the pdf. Consideration of additional higher-order moments (up to order 6) was found insufficient in this regard.

Thus, we conclude from this study that the main challenge with traditional macroscopic models involving moments of the orientation pdf lies more with representation capabilities in highly confined conditions than with a suitable description of the induced orientation kinematics. Use of the averaged fully-confined macroscopic descriptor ˜ a C proposed in this paper is recommended in future theoretical developments. A closed solution cannot be derived, but approximate expressions can be obtained by closing the integral terms. For that purpose, we define These equations require appropriate closure relations for the average higher-order orientation moments ˜ B and ˜ B .

Fig. 1 .

 1 Fig. 1. Hydrodynamic and contact forces acting on a confined rod immersed in a simple shear flow.

Fig. 2 .

 2 Fig. 2. Orientation evolution for H/L = 0 . 8 . Confined (red curve) versus unconfined (broken blue curve) trajectories.

Fig. 3 .

 3 Fig. 3. Orientation evolution for H/L = 0 . 8 (top-left), H/L = 0 . 6 (top-right), H/L = 0 . 4 (bottom-left) and H/L = 0 . 2 (bottom-right).

Fig. 4

 4 Fig.4. Confined Jeffery solution a disc ( t ) versus fully-confined macroscopic orientation tensor a C ( t ) in the case of rods aligned in the same direction. It can be noticed that the planar components of the orientation tensor evolve smoothly whereas the out-of-plane component shows a localized evolution when, after a period of time during which the fibre beads slide at the wall, they detach from it.

Fig. 5 .

 5 Fig.5. Fully-confined Jeffery solution a disc , fc (fibre detachment from the wall is prevented) versus fully-confined macroscopic description a C ( t ) in the case of rods aligned in the same direction.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Isotropic initial orientation distribution for H/L = 0 . 2 and N = 40 0 0 .

Fig. 8 .

 8 Fig. 8. Initial in-plane isotropic orientation distribution related to the fully-confined simulation on z = H, H/L = 0 . 2 , implying p z = 0 . 2 for all fibres.

Fig. 9 .

 9 Fig. 9. Confined Jeffery's orientation tensor a disc ( t ) versus fully-confined macroscopic description a C ( t ) for a population of rods with the fully-confined initial isotropic distribution depicted in Fig. 8 .

Fig. 10 .

 10 Fig.10. Fully-confined Jeffery's orientation tensor versus a C ( t ) for a population of rods with the fully-confined initial isotropic distribution depicted in Fig.8.

Fig. 11 .

 11 Fig. 11.Comparison of confined a disc and unconfined a J , disc solutions computed by using a population of rods large enough to ensure their use as reference solutions, and the orientation tensor a J ( t ) obtained by integrating Eq. (20) using the quadratic and hybrid closure relations, a quad and a hybr respectively, for different initial confined orientation distributions.

Fig. 12 .

 12 Fig. 12. Comparing the orientation development for two different degrees of confinement: H/L = 0 . 2 (top) and H/L = 0 . 3 (bottom): lubrication fully-confined macroscopic ˜ a C (t; H) versus confined Jeffery's a disc ( t ).

1 1 -p 2 z, δ 2 = p z and δ 3 = p z 1 -p 2 z= δ 1 δ 2 .

 122122 ) δ 3 δ 2 b dz • g , (C.2)with p z = z/L and δ 1 =

.(C. 5 )

 5 z) b dz = 0 H 0 ω(z) c dz = 0 H 0 ω(z) B dz = 0 H 0 ω(z) B dz = 0The co-factors involving the delta coefficients in Eqs. (C.1) and (C.2) are noted for the sake of notational simplicity asD 1 = δ 1 , D 2 = δ 3 , D 3 = δ 1 δ 2 , D 4 = δ 22 and D 5 = δ 2 δ 3 . They accept the decompositionD i = ˜ D i + D i , i = 1 , . . . , 5 , (C.6) with H 0 ω(z) D i dz = 0 , ∀ i .Thus, neglecting integrals involving products of variations, as for example H 0 ω(z) D 1 b dz ≈ 0 , and taking into account that tilda variables do not depend on the z -coordinate and thatH 0 ω(z) dz = 1 , Eqs. (C.1) and (C.2) read: ˙ ˜ b = G • ˜ b + ˜ b • G T -2 ˜ D 1 G : ˜ B + (g ⊗ ˜ c + ˜ c ⊗ g ) -

Table 1

 1 Orientation time delay for different degrees of confinement H / L .

	H / L	t/ t J f l
	0 .8	1 .3
	0 .6	1 .9
	0 .4	2 .6
	0 .2	4 .1

  Fully-confined solutions a C ( t ) for different degrees of confinement H / L .

		1
											a xx	(H/L=0.1)
		0.9									a yy	(H/L=0.1)
		0.8									a xx (H/L=0.2) a yy (H/L=0.2)
		0.7									a xx (H/L=0.3)
											a yy	(H/L=0.3)
		0.6									a xx (H/L=0.4)
	a ij	0.5									a yy (H/L=0.4) a xx (H/L=0.5)
		0.4									a yy (H/L=0.5)
		0.3									a xx (H/L=0.6) a yy (H/L=0.6)
		0.2
		0.1
		0	0	1 0	2 0	3 0	4 0	5 0	6 0	7 0	8 0	9 0	1 0 0
								time (s)
	Fig. 13.
													J,disc zz	disappears
													when confinement effects act.

7 , one notices that the orientation a disc xx experiences a slight delay with respect to a J,disc xx , and the peak in a

Acknowledgements

A. Scheuer is a Research Fellow of the "Fonds de la Recherche Scientifique de Belgique" -F.R.S. -FNRS.

Appendix A. Advanced mesoscopic modelling

For representing accurately the orientation distribution at the mesoscopic scale, we could consider two pdf's, ψ J (x , t, p ∈ J ) and ψ C (x , t, p ∈ C) , with

where the normalization condition reads

Now, assuming homogeneous flow for the sake of simplicity and without loss of generality, the dependence of both pdf's on the space coordinates x can be ignored and the Fokker-Planck equation reads:

where, due to symmetry considerations, we only consider the up- and˙ p M represents the velocity on the manifold, defined from

with ˙ p given by Eq. ( 14) .

Here, Q + = ( ˙ p + • t ) ψ J (p ) represents the unconfined rods reaching the manifold C + , t being the unit vector tangent to the unit sphere S, normal to the manifold C + and pointing outward of the allowed region J , with the upstream velocity ˙ p + given by

On the other hand, confined rods leaving the manifold C + are given by Q -, with

The only boundary condition to be prescribed at the boundary of J to ensure conservation of probability reads

where ∂J -denotes the part of the boundary of J through which rods leaving the manifold C come into domain J .

Note that in the Fokker-Planck model (A.3) , domains C and J exchange rods while ensuring conservation of probability. Domain J , due to its 2D nature, exchanges rods through its boundary, whereas C being 1D, the rod exchange appears as a source term in the balance equation (in fact C as previously defined is unbounded).

The numerical treatment of the resulting mesoscopic model is quite delicate because rods leaving the manifold C usually group on the two trajectories in J consisting of the unconfined Jeffery orbits tangent to the manifold C that implies a Dirac delta distribution in J . For this reason, discretization based on the use of continuous approximations remains extremely difficult even when considering two pdf's. The use of a particle-based integration technique constitutes however a plausible route.

Appendix B. Evolution equation for the confined orientation tensor

In this appendix, we address the obtention of Eq. (36) . For that purpose, we consider the confined Jeffery equation

where p is written as

and the gradient of velocity as

The first term of Eq. (B.1) involves

Operating on the second term of Eq. (B.1) , we obtain

Thus, we finally obtain

that can be rewritten as ˙ q = G • q -δ 1 q T • G • q q + δ 2 g -δ 3 q T • g q (B.9)

with δ 1 = 1


Massinissa Merabet

Erian

Miklos Molnar

Generalization of the Minimum Branch Vertices Spanning Tree Problem

Keywords: Spanning Tree, Minimization of Branch Vertices, ILP, Linear Programming, MBVST, k-MBVST, Optical Networks

Given a connected graph G, a vertex v of G is said to be a branch vertex if its degree is strictly greater than 2. The Minimum Branch Vertices Spanning Tree problem (MBVST) consists in nding a spanning tree of G with the minimum number of branch vertices. This problem has been well studied in the literature and has applications specially for routing in optical networks. In this paper we propose a generalization of this problem. We introduce the notion of k-branch vertex which is a vertex with degree strictly greater than k + 2. The parameter k can be seen as the limit of the capacity of optical splitters to divide the light signal. In order to respect as far as possible this limit, we propose to search a spanning tree of G with the minimum number of k-branch vertices (k-MBVST problem). We propose a proof of NP-hardness of this new problem whatever the value of k. We also propose an ILP formulation of the k-MBVST by generalising the MBVST one. Experimental tests on random graphs show that the number of k-branch vertices increases with graph size but decreases with k as well as with the density. They also show that when k ≥ 4, the number of k-branch vertices is close to zero whatever the size and the density of the tested graph.

Introduction

Wavelength-Division Multiplexing (WDM) is an eective technique to exploit the large bandwidth of optical bers to meet the explosive growth of bandwidth demand in the Internet [START_REF] He | Multicasting in wdm networks[END_REF].

Multicast consists in simultaneously transmit information from one source to multiple destinations in a bandwidth ecient way (it duplicates the information only when necessary). From the computational point of view, multicast routing protocols in WDM networks is mainly based on light-trees [START_REF] Sahasrabuddhe | Light trees: optical multicasting for improved performance in wavelength routed networks[END_REF]. It requires the intermediate nodes to have the ability to split the input signal to multiple outputs if needed. A light-splitting switch is needed in the optical device to perform such a task. A node which has the ability to replicating any input signal on any wavelength to any subset of output bers is referred to as a Multicast-Capable (MC) node [START_REF] Malli | Benets of multicasting in all-optical networks[END_REF]. On the other hand, a node which has the ability to tap into the signal and forward it to only one output is called a Multicast-Incapable (MI) node. The light-splitters switchs are rather expensive devices therefore the number of MC nodes should be minimum in the light-tree. L. Gargano et al. expressed in [START_REF] Gargano | Spanning Trees with Bounded Number of Branch Vertices[END_REF] this aspect by introducing The Minimum Branch Vertices Spanning Tree problem (MBVST) which consists in nding a spanning tree of a graph with the minimum number of branch vertices (vertices with degree strictly grater than 2). This NP-hard and no-APX problem [START_REF] Gargano | Spanning Trees with Bounded Number of Branch Vertices[END_REF] is extensively studied in the literature. In [START_REF] Merabet | Minimization of branching in the optical trees with constraints on the degree of nodes[END_REF] that the set of optimal solutions for MBVST and the set of optimal solutions for MDST are disjoint. They also propose two variants of them, taking into account the position of MC nodes in the optical network. In [START_REF] Merabet | Exact solution for branch vertices constrained spanning problems[END_REF], they consider the case where the application do not explicitly impose a sub-graph as solution. A more exible structure called hierarchy is proposed. Hierarchy, which can be seen as a generalization of trees, is dened as a homomorphism of a tree in a graph [START_REF] Molnár | Hierarchies for Constrained Partial Spanning Problems in Graphs[END_REF]. Since minimizing the number of branch vertices in a hierarchy does not make sense, they propose to search the minimum cost spanning hierarchy such that the number of branch vertices is less than or equal to an integer R.

The light power in optical networks should be controlled because of the power loss. Indeed, if a light signal is splitted into k copies, the signal power of one copy will be reduced with, at least, a factor of 1/k of the original signal power [START_REF] Maher | Power-Ecient Design of Multicast Wavelength-Routed Networks[END_REF]. The more k is high, the more informations are lost. If k is too high then the information can not be read by the destinations. Therefore k can be seen as a tolerance parameter. Let a k-branch vertex be a vertex with degree strictly grater than k + 2, it is useful to look for a light-tree in the WDM network with the minimum number of k-branch vertices, where k is xed as the tolerance parameter. Notice that the MBVST is a special case with k = 0. If the light-tree contains k-branch vertices, an optical amplier should be installed near each one to guaranty the eciency of the broadcast/multicast.

In this paper we introduce the k-MBVST problem which, given a positive integer k and a graph G = (V, E), aims to nd a minimum spanning tree of G with minimum number of k-branch vertices. I proof of the NP-hardness of the k-MBVST whatever k is given in Section 2. An ILP formulation based on single commodity ow is given in Section 3, and the experimentations on random graphs are done in Section 4.

Problem formulation and NP-hardness

The MBVST problem is dened in [START_REF] Gargano | Spanning Trees with Bounded Number of Branch Vertices[END_REF] as follows :

Denition 2.1 Let G = (V, E) be a graph. The MBVST problem consists in nding a spanning tree T of G such that the number of branch vertices in T is minimum.

Let us dene the k-branch vertex : Denition 2.2 A k-branch vertex is a vertex with degree strictly grater than

k + 2.
The k-MBVST problem which can be seen as the generalization of the MBVST problem is dened as follows :

Denition 2.3 Let G = (V, E) be a graph and k be a positive integer. The k-MBVST problem consists in nding a spanning tree T of G such that the number of k-branch vertices in T is minimum. We denote by s k (G) the smallest number of k-branch vertices in any spanning tree of G.

In a Hamiltonian graph, nding a 0-MBVST is equivalent to nding a Hamiltonian path in G. Thus, The k-MBVST is NP-complete in this case. Furthermore, the classical MBVST problem is NP-complete [START_REF] Gargano | Spanning Trees with Bounded Number of Branch Vertices[END_REF] and it is a particular case of the k-MBVST problem. Therefore k-MBVST is at least as dicult as the MBVST. We proof in the following that the k-MBVST problem still hard even when k is distant from zero.

Theorem 2.4 Let r be a xed non-negative integer. It is NP-complete to decide whether a given graph G satises s k (G) ≤ r whatever the value of k.

Proof

• For r = 0 : Let G = (V, E) be a given graph. Construct a new graph G by linking k leafs to each vertex v ∈ V . Decide whatever G contains a spanning tree with no k-branch vertex is equivalent to decide if G is Hamiltonian. • For r ≥ 1 : Let G = (V, E) be a given graph. We construct a graph G by duplicating r • k k times the graph G and adding a complete graph K r .

We choose an arbitrary vertex v ∈ V and we link every vertex of the graph K to r distinct duplications of G from their vertex v. We link k leafs to each vertex of each duplication of G. In any spanning tree of G , the k vertices of the complete graph are necessary a k-branch vertices. Thus, the graph G contains a spanning tree with s k (G ) = r if and only if G admits a Hamiltonian path starting from v. In this section, we propose an ILP formulation of the k-MBVST problem based on the single ow formulation proposed in [START_REF] Cerulli | Bounded-degree spanning tree problems: models and new algorithms[END_REF]. In order to dene a spanning tree T of G, we can send from a source vertex s ∈ V one ow unit to every other vertices v ∈ V \ {s}. The ow is directed. Thus, the initial graph has to be transformed in a symmetric oriented graph G d = (V, E d ) where for each edge {u, v} in E corresponds two arcs (u, v) and (v, u) in E d . For each arc e = (u, v) ∈ E d , we dene a variable f (u,v) representing the ow going from u to v. We have a binary decision variable x e . It is equal to 1 if f (u,v) or f (v,u) carry out a non-zero ow, and zero otherwise. Finally, for each v ∈ V , we have a decision variable y v that is equal to 1 if v is a k-branch vertex, and 0 otherwise. In the following the linear program is presented.

v (a) Graph G v v v v v v v v v v v v v v v v v v (b) Graph G
The objective of our problem is to minimize the number of k-branch vertices belonging to the spanning tree of G. Hence the general objective function can be expressed as follows:

M inimize : v∈V y v (1)
This objective function is subject to a set of constraints.

Spanning tree constraints:

(u,v)∈E d x (u,v) = 1 ∀v ∈ V \ {s}
(2)

(u,v)∈E d x (u,v) = n -1 (3)
Since a vertex with more than one parent creates a cycle, the constraint (2) ensure that each vertex except the source have one and only one predecessor.

The number of edges of any tree must be equal to the number of its vertices minus one. The constraints (3) makes sure that exactly n -1 arcs are selected in the solution. These two constraints are necessary but not sucient to have a tree since the connectivity must be ensured. In this purpose, ow based constraints are added.

Connectivity constraints:

(s,v)∈E d f (s,v) - (v,s)∈E d f (v,s) = |V | -1 (4) (v,u)∈E d f (v,u) - (u,v)∈E d f (u,v) = -1 ∀v ∈ V \ {s}
(5)

x (u,v) ≤ f (u,v) ≤ (|V | -1) • x (u,v) ∀(u, v) ∈ E d
(6) Constraints (5) ensures that each vertex except the source "consumes" one and only one unit of ow. This constraint also guarantees that each vertex is reachable from the source s. Constraint (4) ensure that the ow emitted by the source is equal to |V | -1. Constraints (6) allows each arc to carry non-zero ow if and only if it is used in the output graph. The value of this ow should not exceed the total ow emitted by the source.

Degree constraints:

(v,u)∈E d x (v,u) + (u,v)∈E d x (u,v) -k -2 ≤ d(v) • y v ∀v ∈ V (7)
Constraints (7) impose vertex v to be a k-branch vertex if its degree is strictly greater than k + 2 in the tree. Note that the value of y v is unconstrained if d(v) ≤ k + 2, however in this case it will be set to zero by the objective function.

Experimentation

In this section, we describe the computational results that we obtained by applying the proposed single commodity ow formulation (SC) for the k-MBVST to a set of instances generated according to parameters originally proposed in [START_REF] Carrabs | Lower and upper bounds for the spanning tree with minimum branch vertices[END_REF]. In order to obtain a signicant number of branch vertices, these instances are considered as spars. We consider 9 dierent values for the number of vertices of random graph: |V | = {50, 100, 200, 300, 400, 500, 600, 700, 800}.

The number of edges is generated according to the following formula:

(|V | -1) + i × 1.5 × |V | with i ∈ {1, 2, 3}.
For each value of the parameter k ∈ {0, 1, 2, 3, 4, 5}, we randomly generated 30 instances for each choice of |V | and i. The SC formulation has been coded in C. The IBM ILOG CPLEX 12 solver was used to solve the mathematical formulations, considering a time limit of 1 hour for each instance. All tests have been executed on an Intel i7 6820HQ 2.7Ghz (with 8 Cores) workstation with 16 gigabytes of RAM.

The numerical results are presented in Table 1. Trivially, the computational time increases both with the size of the graph and with de density. Indeed, greater and denser is the instance, higher is the number of decision variables in the ILP. Furthermore, greater is the instance, less the constraint 6 is tight in average. The number of branch vertices decrease in a more dense instances (i = 2 and i = 3) because more an instance is dense and more farther it is from the Hamiltonicity form. Obviously, the number of branch vertices mechanically increase with the instance size. This is amplied by the fact that our density formula makes the density decrease with the instance size.

Figure 3 shows the number of k-branch vertices regarding the variation of the parameter k and |V | for each value of i. The number of k-branch vertices mechanically increases with |V | but decreases with k as well as with i. When k ≥ 4, the number of k-branch vertices is close to zero whatever the value of |V | and the value of i. This armation can not be conrmed for the unsolved instances but even when the percentage of solved instances is equal to 100% this armation still true.

k=0 Instances i = 1 i = 2 i = 3 |V | Sol Time # Inst

Conclusion

In this paper, we propose a generalization of the well known MBVST problem by introducing the notion of k-branch vertex, which is a vertex with degree strictly greater than k + 2. Our new parametrized problem (k-MBVST) aims to nd a spanning tree of G with the minimum number of k-branch vertices.

Let r be a non-negative integer, we proved that it is N P -complete to decide whatever a graph can be spans by a tree with at most r k-branch vertices, whatever the value of k. We also proposed an ILP based on a single ow formulation. Tests on spars random graphs allowed us to evaluate the number of k-branch vertices in the optimal solution and the running time regarding the value of the parameters k, the graph size, and the graph density. The results show that the number of k-branch vertices increases with graph size but decreases with k as well as with the density. They also show that when k ≥ 4, the number of k-branch vertices is close to zero whatever the size and the density of the tested graph.

  [CGI09], Cerrulli et al. give the rst ILP formulation of this problem based on a single commodity ow to guaranty the connectivity. In [CCGG13] F.Carabbs et al. give two more ILP formulations based respectively on Multi Commodity Flow and Miller-Tucker-Zemlin formulation. They also provide both lower and upper bound for the MBVST using Lagrangian relaxation. In [Mar15], A.Marin presents a branch-and-cut algorithm based on an enforced Integer Programming formulation for the MBVST problem. In [CCR14], C.Cerrone et al. present a unied memetic algorithm for the MBVST, the problem of minimize the degree sum of branch vertices (MDST), and the well known Minimum Leaves Problem. M.Merabet et al. prove in

Figure 1 .

 1 Figure 1. Reduction from Hamiltonien problem to 0-MBVST (k = 5)

Figure 2 .

 2 Figure 2. Construction of the graph G for k = 1 and r = 6.

Figure 3 .

 3 Figure 3. Variation of s k (G) regarding |V |, k and i.

Table 1 Solution

 1

					Sol	Time	# Inst	Sol	Time	# Inst
	50	8.00	0.22	30	3.80	0.60	30	2.50	0.62	30
	100	18.00	0.39	30	12.30	1.11	30	7.60	5.44	30
	200	38.40	1.14	30	28.90	5.21	30	20.70	18.31	30
	300	59.90	1.91	30	49.40	9.73	30	37.90	42.76	30
	400	80.60	4.06	30	67.40	15.15	30	53.80	50.98	30
	500	102.60	4.32	30	86.60	19.69	30	74.10	446.47	30
	600	130.50	5.93	30	107.70	21.81	30	91.22	289.46	29
	700	149.80	6.68	30	131.70	40.04	30	113.44	275.39	29
	800	175.20	8.73	30	149.80	43.79	30	131.00	188.82	30
					k=1
	50	1.50	0.10	30	0.00	0.21	30	0.00	0.25	30
	100	4.60	0.31	30	1.00	0.46	30	0.50	0.82	30
	200	12.60	0.58	30	5.30	1.80	30	2.67	3.73	27
	300	20.20	1.72	30	10.70	8.32	30	4.88	9.35	24
	400	28.10	3.21	30	17.70	12.19	30	9.50	21.20	25
	500	41.00	4.78	30	24.67	16.92	27	14.00	42.85	27
	600	49.00	6.16	30	31.89	22.37	27	20.43	108.02	21
	700	61.20	8.11	30	38.90	23.84	30	31.00	77.04	6
	800	68.70	8.99	30	46.60	40.74	30	34.60	150.90	12
					k=2
	50	0.30	0.07	30	0.00	0.12	30	0.00	0.17	30
	100	0.70	0.16	30	0.10	0.31	30	0.00	0.49	30
	200	3.00	0.49	30	0.50	1.45	30	0.33	1.68	27
	300	4.90	0.94	30	1.80	3.56	30	0.71	6.19	21
	400	9.80	2.09	30	2.67	7.36	27	1.22	15.07	27
	500	12.00	2.29	30	4.63	12.04	24	2.50	25.04	30
	600	16.20	4.23	30	6.33	24.68	18	3.25	32.69	12
	700	20.00	5.56	30	9.75	19.19	12	5.00	118.11	6
	800	24.70	5.62	30	12.00	23.90	9	4.33	57.01	9
					k=3
	50	0.10	0.07	30	0.00	0.10	30	0.00	0.13	30
	100	0.10	0.15	30	0.10	0.17	30	0.00	0.27	30
	200	0.60	0.32	30	0.10	0.93	30	0.00	1.56	30
	300	1.10	0.79	30	0.50	1.63	30	0.10	3.98	30
	400	1.40	1.32	30	1.00	4.42	24	0.40	10.80	30
	500	3.50	1.62	30	1.38	9.50	24	0.13	22.14	24
	600	4.56	3.00	27	2.00	5.59	7	1.00	30.55	21
	700	6.20	4.18	30	1.40	53.08	15	0.33	69.49	9
	800	7.30	3.04	30	2.14	21.57	20	2.50	57.01	6
					k=4
	50	0.00	0.05	30	0.00	0.09	30	0.00	0.12	30
	100	0.00	0.08	30	0.00	0.13	30	0.00	0.23	30
	200	0.00	0.18	30	0.00	0.54	30	0.00	0.35	30
	300	0.50	0.32	30	0.14	0.98	21	0.00	1.88	27
	400	0.60	0.54	30	0.00	09.12	21	0.00	2.41	24
	500	0.40	0.51	30	0.13	5.00	24	0.00	24.68	18
	600	1.60	1.15	30	0.20	8.65	15	0.25	50.22	12
	700	1.30	2.63	30	0.67	15.77	9	0.00	86.80	9
	800	1.80	3.59	30	0.33	16.36	18	0.40	14.20	15
					k=5
	50	0.00	0.05	30	0.00	0.08	30	0.00	0.10	30
	100	0.00	0.09	30	0.00	0.15	30	0.00	0.14	30
	200	0.00	0.11	30	0.00	0.42	30	0.00	0.50	30
	300	0.00	0.20	30	0.00	1.12	24	0.00	1.08	30
	400	0.30	0.34	30	0.00	21.68	24	0.00	4.76	24
	500	0.10	0.29	30	0.00	80.04	21	0.00	45.24	18
	600	0.30	0.64	30	0.00	6.34	15	0.00	40.58	15
	700	0.00	1.83	21	0.00	16.95	8	0.00	42.75	8
	800	0.11	2.39	27	0.00	10.39	21	0.00	49.27	6

value, running time and number of solved instances regarding |V |, k and i.

Email: mmerabet@ntu.edu.sg

Email: molnar@lirmm.fr