On the multiscale description of dilute suspensions of non-Brownian rigid clusters composed of rods

The motion of an ellipsoidal particle immersed in a ﬂow of a Newtonian ﬂuid was obtained in the pio-neering work of Jeffery in 1922. Suspensions of industrial interest usually involve particles with a variety of shapes. Moreover, suspensions composed of rods (a limit case of an ellipsoid) aggregate, leading to clusters with particular shapes that exhibit, when immersed in a ﬂow, an almost rigid motion. In this work, we propose a framework for describing dilute suspensions of rigid particles and derive an expression for calculating the motion of rigid clusters of general shape immersed in a ﬂow of a Newtonian ﬂuid. We show that the cluster’s rotary velocity only depends on a symmetric tensor c with unit trace that can be considered as the appropriate conformation tensor for describing cluster kinematics.


Introduction
Suspensions involving particles can be described at the microscopic scale by tracking the motion of each one of the particles involved in the system.This approach is based on three main elements: (i) the knowledge of the equation governing the particle motion in the fluid flow; (ii) the introduction of the particle effects on the flow kinematics if coupled simulations are envisaged; and (iii) the availability of computational resources for tracking efficiently millions of particles.
In dilute suspensions, the motion of ellipsoidal particles can be accurately described by using Jeffery's equation [29].When the concentration becomes large enough, interactions cannot be neglected any longer and the calculation becomes more complex from the computational point of view.At this scale, currently available simulations remain quite far from the scenarios of industrial interest.
For circumventing the difficulties related to simulations at the microscopic scale, these being more computational than conceptual, coarser models were introduced.The interested reader can refer to [31] and the references therein for a review on multiscale approaches in the context of computational rheology.
Mesoscopic kinetic theory models result from coarsening microscopic descriptions.In kinetic theory models the individuality of the particles is lost in favor of a statistical description that substitutes the entities by a series of conformation coordinates [21,7].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 the probability distribution function -pdfwhose conservation balance results in the so-called Fokker-Planck equation governing its evolution.Fokker-Planck equations involve the flow induced conformation evolution.In the case of a suspension of rods, the flow induced conformation (orientation) evolution is given, as indicated above, by Jeffery's equation.Since the pdf depends on the physical coordinates (space and time) and a series of conformational coordinates, the associated Fokker-Planck equation is multidimensional.Standard mesh-based discretization techniques fail when addressing multidimensional models.This issue is known as the curse of dimensionality and it justifies the few number of existing works addressing the solution of kinetic theory models within the Fokker-Planck framework.
For circumventing the curse of dimensionality at the mesoscopic scale, several techniques based on the use of particles were proposed and widely employed.Here the particles are not real particles, but rather should viewed as computational particles that allow one to describe the main suspension features (rheology, properties related to the particles conformation, etc.).Despite the fact of considering a discrete description, the level of detail in the description and the richness of the physics are exactly the same that the ones associated with the use of Fokker-Planck descriptions, and obviously the solutions computed by using both descriptions are in the limit of convergence exactly the same.
The use of the continuous description based on the solution of the Fokker-Planck equation remains challenging because of the high dimensionality that it involves.On the other hand, when employing its discrete counterpart, the main difficulty is related to the extremely large number of particles to be considered.This number depends on the model output of interest.When only the moments of the distribution are concerned, a moderate number of particles is enough.However, when one is interested in the pdf itself, the number of computational particles could become extremely large.
Solution procedures based on the use of particles at the mesoscopic scale have been extensively employed by many authors [38,11,42,14,19,20].On the other hand, there are few works focusing on the solution of Fokker-Planck equations by using standard discretization techniques [33,12].We proposed some years ago a new solution technique called Proper Generalized Decomposition based on the use of separated representations in order to ensure that the complexity scales linearly with the model dimensionality [4,5].This technique consists in expressing the unknown field as a finite sum of functional products, i.e. expressing a generic multidimensional function uðx 1 ; . . .; x d Þ as: The interested reader can refer to [37,15,17] and the references therein for a deep analysis of this technique and its applications in computational rheology.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 is unpredictable [30,13].Alternatively, macroscopic equations are carefully postulated, within a top-down approach, in order to guarantee the model objectivity and its thermodynamical admissibility.
In the case of dilute suspensions involving ellipsoidal particles, the three scales have been extensively considered starting from the pioneer work by Jeffery in 1922 [29] (see the review by Petrie [39,24,6,23,[25][26][27]). The multiscale modeling involves in this case nine main conceptual bricks [16]: 1. Particle conformation.Choice of the coordinate describing the particle conformation.For ellipsoidal particles, it is usual to consider the unit vector p aligned along the ellipsoid axis.2. Particle conformation evolution.Derivation of the equation governing the time evolution _ p of the particle conformation.For an ellipsoid suspended in a Newtonian viscous fluid, whose kinematics is described by the gradient of velocity rv, one obtains Jeffery's equation [29]: where X ¼ 1 2 rv À ðrvÞ T and D ¼ 1 2 rv þ ðrvÞ T are respectively the skew-symmetric and symmetric components of the gradient of velocity tensor, and k is a coefficient that depends on the ellipsoid aspect ratio r (ratio of length 2L to diameter . In view of the symmetry of the strain rate tensor D, the equality rv : ðp pÞ ¼ D : ðp pÞ applies.
Three special cases can be considered: Slender body theories result for infinite aspect ratio ellipsoids, implying r ! 1 and then k ! 1.Thus Jeffery's equation reduces to: When the ellipsoid length 2L coincides with the diameter of its circular cross section / the ellipsoid reduces to a sphere and in that case r ¼ 1 and k ¼ 0, implying that the sphere rotates according to the fluid flow vorticity X, that is Finally, when the ellipsoid length reduces to zero with respect to its cross section diameter, the ellipsoid becomes a flat and thin disk.In that case r !0 and k !À1.The disk kinematics is then given by where p is the unit vector normal to the disk surface.3. Particle contribution to the stress.4. Population description.Description of a population of N ellipsoidal particles located at position x.Note that these particles have a computational nature and do not correspond to the ones that would have been concerned by a purely microscopic description.Thus, the mesoscopic description can be done from the specification of the conformation of each one, that is p i ; i ¼ 1; . . .; N (discrete description); or from the introduction of a probability distribution function -pdfwðx; t; pÞ that gives the fraction of particles that at position x and time t are aligned along the direction p (continuous description).In the case of homogeneous distributions, the space dependence can be ignored and then wðt; pÞ. 5. Description of the population evolution.The microstructure evolution can be described from the time evolution of each particle, i.e. from the integration of _ p i according to Jeffery's equation for each particle i ¼ 1; . . .; N , or from the integration of the Fokker-Planck equation governing the evolution of the pdf wðx; t; pÞ.This equation expresses conservation of probability.It reads: where r x and r p refer respectively to the gradients in spatial and conformational coordinates.If the vector p is defined from the ellipsoid center of gravity G and inertia effects are neglected, _ x ¼ v as it can be shown that the particle center of gravity G moves with the fluid velocity.We come back to this point later.6. Contribution of the particle population to the stress.7. Microstructural macroscopic description.Microstructural macroscopic descriptions are usually based on the use of some pdf moments [2].In view of the symmetry of the pdf (wðpÞ ¼ wðÀpÞ) in the case of rigid rods addressed here, oddorder moments vanish, and the first non-vanishing moment is the so-called second-order orientation tensor defined as: where S refers to the surface of the unit sphere (unit circle in 2D).
The next non-zero moment is the fourth-order orientation tensor A defined as: 8. Microstructural macroscopic evolution.Finally, the microstructural macroscopic evolution can be described from the time derivative of the different moments.One could expect that by taking the time derivative of Eq. ( 7) and substituting the time derivative of p by the expression given by Jeffery's Eq. ( 2), one could obtain an equation for _ a. Indeed, it is the case as we illustrate later, but unfortunately the evolution of tensor a does not only depend on a itself, but also on A. The evolution equation for A itself depends on the sixthorder orientation tensor and so on.The model must be closed and for that purpose one must invoke approximate closure relations between higher and lower-order moments, for example A ¼ AðaÞ.The pertinence and impact of such closures motivated many works [3,22,32,40].9. Moment based stress.
The above conceptual bricks can be extended to Brownian particles, to semi-dilute or semi-concentrated regimes, and even to concentrated and highly-concentrated regimes involving aggregates.In this work, we address dilute suspensions of non-Brownian rigid clusters composed of rods.
First, we propose an original methodology for the above framework in the case of dilute suspensions of non-Brownian rods (ellipsoids with infinite aspect ratio).In Section 3, we extend that methodology to rigid clusters composed of rods of arbitrary shape.In Section 4, the resulting models are validated in the case of ellipsoids, spheres, disks and cylinders by comparison with existing analytical or numerical solutions.Finally, we conclude in Section 5 and give an overview of perspectives on the above multiscale approach, most of them being work in progress within our group.
Remark 1.In this paper we consider the following tensor products, where Einstein's summation convention is assumed:

Multiscale modeling of dilute suspensions of non-Brownian rods
We consider a suspending medium consisting of a Newtonian fluid of viscosity g 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 that is defined by the velocity field vðx; tÞ, with x 2 X 2 R d ; d ¼ 2; 3.In what follows we consider in detail the different conceptual bricks previously introduced.

Particle conformation
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 rods axis.Inertial effects are neglected in the sequel.

Particle conformation evolution
The orientation evolution of an ellipsoidal particle is described by Jeffery's Eq. ( 2), that in the case of rods with infinite aspect ratio reduces to Eq. ( 3).
Eq. ( 3) can be derived in a very simple way illustrated in Fig. 1.We consider a system 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 þ rv Á pL and the second one by v G þ _ pL.Thus, the force FðpLÞ reads: where n is the friction coefficient, v 0 the fluid velocity at the rod center of gravity, and v G the velocity of the center of gravity.
If F applies on the bead at pL, then the force on the opposite bead at ÀpL reads By adding Eqs. ( 9) and ( 10), and neglecting inertial effects, we obtain the force balance that is, the rod center of gravity is moving with the fluid velocity.For simplicity of notation, we shall write F ¼ FðpLÞ and FðÀpLÞ ¼ ÀF.
As the resulting torque must also vanish, the only possibility is that force F acts along p, that is F ¼ kp, with k 2 R. Thus, we can write Premultiplying Eq. ( 12) by p and taking into account that p Á p ¼ 1 and consequently p Á _ p ¼ 0, we have We have thus obtained Jeffery's Eq. ( 3), Fig. 1.Hydrodynamic forces applied on a rod immersed in a Newtonian fluid.
Remark 2. As the factor nL appears in both sides of Eq. ( 14), the rod kinematics does not contain size effects.

Particle contribution to the stress
The forces applied at the rod ends pL and ÀpL are respectively kp and Àkp, i.e. directed along the rod and in equilibrium by construction.
With k given by Eq. ( 13), we have By applying Kramers' formula, the corresponding contribution to the stress is given by which can be rewritten as Remark 3. From Eq. ( 18), one could infer that rheology contains size effects in view of the presence of factor L 2 in the particle contribution to the stress.To prove that size effects also disappears when addressing rheology, we consider the traction T applied on the rod end.Considering Newtonian behavior, we have By writing the objective stress related to the rod rotation 2gð _ p À X Á pÞ and enforcing that the resulting force aligns along the direction p, we obtain again Jeffery's equation.In order to obtain the same expression for the stress, it suffices to consider n / g L 2 .Thus, we can conclude that size effects are absent in Eq. ( 18).

Population description
As previously mentioned, there are two natural descriptions of a population of rods.
The first one consists in specifying each rod orientation by considering the unit vector aligned along its axis, that is, by considering p i ; i ¼ 1; . . .; N .As discussed in the next section, the main drawback of this approach lies in the necessity of tracking the evolution of each ''computational'' rod by solving the corresponding Jeffery equation, and even if conceptually there is no major difficulty, the computing cost could be excessive in most practical applications.The second approach is the introduction of the pdf wðx; t; pÞ that gives the fraction of rods that a position x and time t are oriented along direction p.
Despite the fact that both mesoscopic models involve the same physics and richness of description, the main advantage of the second one is the manipulation of a scalar continuous function instead of the discrete description involved in the first approach.The price to be paid when using the description based on the use of the pdf is its inherent multidimensionality, because in that framework the pdf depends on the standard space and time coordinates and also on the conformation coordinates that the microstructural description involves, i.e. p in the present case.

Description of the population evolution
When the population is described from the individuals composing it, whose conformation is given by vectors p i ; i ¼ 1; . . .; N , the evolution of each one is given by Jeffery's equation: The alternative description consists of using the pdf wðx; t; pÞ that satisfies the normalization condition: Z S wðx; t; pÞ dp ¼ 1; 8x; 8t: where, for inertialess rods, _ x ¼ vðx; tÞ, and the rod rotary velocity is given by Jeffery's equation: The price to pay is the increase of the model dimensionality, as the orientation distribution is defined in a 6-dimensional domain, i.e. w : 2.6.Contribution of the particle population to the stress Again, we consider the two alternative descriptions: When the population is described in a discrete manner by means of the p i vectors, the contribution of rods to the suspension stress is calculated by adding their individual effects given by Eq. ( 18), that is: where N ðx; tÞ refers to the number of computational rods located at time t in the neighborhood of x.When the population is described from the pdf, the sum in Eq. ( 23) is replaced by an integral in the conformation space S: Here, the particle number N p accounts for the particle concentration and the viscosity is used instead of the friction coefficient to be consistent with the usual notation.
In terms of the fourth-order orientation tensor (8), we obtain: which in view of the symmetry of A can be rewritten as

Macroscopic description
As just discussed, discrete descriptions are computationally expensive because of the large number of rods that must be considered in order to derive accurate-enough model outputs.On the other hand, Fokker-Planck descriptions are rarely considered in view of the curse of dimensionality that introduction of conformation coordinates implies.Thus, standard mesh-based discretization techniques, as finite differences, finite elements or finite volumes, fail when addressing models defined in high-dimensional spaces.
For these reasons, mesoscopic models were coarsened to derive macroscopic models defined in standard physical domains, involving space and time.At the macroscopic scale, the orientation distribution function is substituted by its moments for describing the microstructure.Usually, macroscopic descriptions of rod suspensions are based on the use of the first two non-zero moments, i.e.
the second and the fourth-order moments defined in Eqs. ( 7) and ( 8) respectively.

Microstructural macroscopic evolution
The microstructural evolution described at the macroscopic scale considers the time evolution of the pdf moments.The time evolution of the second-order orientation tensor is given by: As discussed previously, this equation involves the fourth-order moment A. The time derivative of the fourth-order moment, using the same rationale, involves the sixth-order moment A, and so on.Thus, an approximate closure relation is needed in order to express the fourth-order moment A as a function of the lowerorder moment a. Different closure relations have been introduced and widely used [3,22,40,32].With the quadratic closure relation (that is only exact when all rods are locally aligned in the same direction), the fourth-order moment is approximated as follows: This gives _ a % rv Á a þ a Á ðrvÞ T À 2 ðrv : aÞ a; We obtained previously the expression of the rod population contribution to the stress: s ¼ 2gN p ðA : rvÞ; ð31Þ which involves the fourth-order moment A. There is no closure issues when A is calculated from the pdf w by using (8).When one proceeds at the macroscopic scale, however, wherein the pdf is not available, a closure relation must be considered for either writing A from the knowledge of a, itself being calculated by integrating (27) with an appropriate closure relation (e.g.Eq. ( 29) when considering the quadratic closure), or calculating A by solving the equation that governs its time evolution in which, as just commented, the sixth-order moment appears requiring again an appropriate closure.
The first route is the simplest one and the most used in practice.It leads to s ¼ 2gN p ðA cr ðaÞ : rvÞ; ð32Þ where the superscript cr refers to the use of an appropriate closure relation.
With the quadratic closure, the stress reads:

Extension to non-Brownian rigid clusters of general shape and composed of rods
After having considered rods immersed in a flow, we extend the above approach to more complex configurations.Suspensions of industrial interest composed of rods (a limit case of an ellipsoid) aggregate, leading to clusters with particular shapes that exhibit, when immersed in a flow, an almost rigid motion.These situations are currently encountered when considering carbon nanotubes suspensions as discussed in [34][35][36].
We represent rigid clusters of general shape in both a discrete manner, assuming they are all composed of N=2 rods (N being an even integer) involving N beads, and in a continuous manner from the continuous pdf describing the configuration of those rod beads.In a rigid cluster, there is no relative motion between the rods composing it.

Discrete description
First, we consider a 3D rigid cluster consisting of N=2 rods R j of length Q j .We assume that each rod R j contains two beads at its ends on which hydrodynamic forces apply.Thus, the cluster contains N beads B i ; i ¼ 1; . . .; N. The location of each bead B i with respect to the cluster center of gravity G is represented by L i p i , where p i is the unit vector pointing from G to B i .The cluster is sketched in Fig. 2.
Brownian effects are neglected and then only flow-induced hydrodynamic forces must be considered (rod-rod hydrodynamic interactions are not considered in this work).Forces F i apply on each bead located at position L i p i (Fig. 2) and are proportional to the difference of velocity between the one of the flow unperturbed by the presence of the cluster at the bead location and the one of the bead itself: where n is the friction coefficient, v the flow velocity field, v 0 the fluid velocity at position G and v G the velocity of the center of gravity.
By adding all the forces we obtain Both sums in Eq. ( 35) vanish, the first one as a direct consequence of the definition of the center of gravity, and the second because the cluster is assumed rigid.Thus, Eq. ( 35) becomes implying that the cluster center of gravity is moving with the fluid velocity at that position.The torque created by forces applied on bead i is given by Neglecting inertial effects, the resulting torque for the whole cluster must vanish: Taking Eqs. ( 37) and ( 34) into account, we have If we define the cluster angular velocity x such that torque equilibrium reads Now, by using the vector triple product relationship a Â ðb Â cÞ ¼ b ða Á cÞ À c ða Á bÞ, the right-hand side reads which, taking into account the normality of vectors p i and the fact with I the unit matrix.Now, the left-hand side of Eq. ( 41) can be rewritten by using the third-order Levi-Civita permutation tensor such that ðu Â vÞ ¼ : ðv uÞ.We obtain which finally yields Expressed in a more compact form, we have

Continuous description
The continuous description considers the pdf !ðp; LÞ giving the fraction of rod beads located at position pL, that can be expressed as !ðp;LÞ ¼ wðp; LÞ CðLÞ; In particular Eq. ( 46) becomes: and due to the fact that trðc L Þ ¼ 1, we can infer that b is in fact the trace of c, that is b ¼ trð cÞ.Thus, the particle kinematics finally results: or dividing by trð cÞ, where the conformation tensor c is given by Remark 4. In the 2D case, c Á x ¼ 0 and then Þ .Moreover, if all beads are located at the same distance L, we recover the expression introduced in [1].
An extremely important consequence of this analysis is that rigid clusters composed of rods having the same conformation tensor c rotate at the same angular velocity.
As the conformation tensor c is symmetric and positive definite, it has real eigenvalues and eigenvectors.In 3D, the three mutually perpendicular eigenvectors will be denoted by u 1 ; u 2 and u 3 , with the associated eigenvalues s 1 ; s 2 and s 3 respectively.Imagine a rigid cluster composed of three rods oriented in directions u 1 ; u 2 and u 3 with respective lengths ffiffiffiffiffi The conformation tensor of such a three-rod cluster coincides with c and then both tensors have the same rotary velocity.
If two eigenvalues are identical, e.g.s 2 ¼ s 3 , it can be shown using the same rationale as in Appendix B that the three-rod cluster rotary velocity coincides with the one associated to an ellipsoid of aspect ratio r ¼ ffiffiffiffi given by Jeffery's Eq. ( 2).

Ellipsoidal clusters
We prove in Appendix B that the motion of a 2D cluster composed of two perpendicular rods of lengths L 1 and L 2 coincides with the Jeffery motion of an ellipse of aspect ratio L 1 =L 2 .This results can be easily generalized to 3D.In that case the motion of a cluster composed of three mutually perpendicular rods of lengths L 1 and L 2 for two of them, coincides with the motion predicted by Jeffery for an ellipsoid [29].This equivalence can be verified for any aspect ratio, for r ¼ L 1 L 2 ranging from r ¼ 0 to 1 : r ! 1 (rods), r > 1 (prolates), r ¼ 1 (sphere), 0 < r < 1 (oblates) and r !0 (disks).We would like to emphasize that this is not an approximate result, it is a formal equivalence.
In the case of a non-axisymmetric ellipsoid, the motion of a cluster composed of three mutually perpendicular rods of lengths L 1 ; L 2 and L 3 coincides with the one predicted by Hinch and Leal [28] for a non-axisymmetric ellipsoid having the same axis lengths.
As shown previously, the cluster rotary velocity only depends on the conformation tensor, which is symmetric and with unit trace.The cluster can be reduced to two or three rods, in 2D and 3D respectively, whose lengths and orientations are given by the eigenvalues and eigenvectors of the conformation tensor c.
This equivalence between the motion of an ellipsoid and of an equivalent three-rod cluster may appear strange to some and natural to others.Now the question is: what happens if starting from the three-rod cluster (representing its equivalent ellipsoid) we add new rods to it?
Imagine that the added rods have a length such that their terminal beads are located on the surface of the associated ellipsoid.One could expect that in such circumstance the cluster motion remains unchanged.However, it is easy to understand that the addition of the new rods will modify the conformation tensor c and then its kinematics.The motion of the new resulting cluster will be the same as that of an equivalent ellipsoid (defined from the eigenvalues and eigenvectors of c) but which is no longer the imaginary-one on which the beads of the cluster are located.
There is however a particular case that is quite curious.Imagine a (non necessarily axisymmetric) ellipsoid filled by a uniform distribution of beads.This situation is not very physical because in such circumstances the flow will be extremely perturbed at the bead positions located inside due to the screening of more external layers of beads.Now, the motion associated with the resulting conformation tensor c coincides with Jeffery's motion of that ellipsoid.Again, this is not an approximate result: there exists a formal equivalence for any ellipsoid.

Cylindrical clusters
In this section, we consider the case of more general clusters shapes.
First, we consider an hypothetical cylinder of diameter D and height H, with aspect ratio r ¼ H D .Many authors tried to define equivalent ellipsoid aspect ratios for cylinders [10,9,18,8,25,41].The interested reader can refer to Section 2.2 in [39].Different equivalent ellipsoid aspect ratios were given for shape-ended and blunt-ended cylinders.An interesting recent work was carried out by Zhang et al. [43].The authors determined by using 3D finite element simulations the equivalent ellipsoid aspect ratio r e for cylinders with different aspect ratio r.They proposed the relationship: We consider now different rigid clusters composed of rods with their center of gravity located at G ¼ ð0; 0; 0Þ, and their terminal beads located on the surface of an hypothetical cylinder.Thus, we consider N Á M beads B i;j , with coordinates . . .; N and j ¼ 1; . . .; M), N and h z ¼ H MÀ1 .We study different configurations by selecting different values of N and M, computing the associated conformation tensors c from c ¼ and from them the equivalent ellipsoid aspect ratios r e that will be compared with the predictions given by Eq. ( 57).Results are summarized in Table 1.Then, we compare motions predicted by Jeffery's model for ellipsoids having equivalent aspect ratios given by the Zhang et al. model (57) and from the cluster model (46).Table 1 shows that N does not influence the results.This fact can be easily understood because the circular shape of the ellipsoid cross section is perfectly described with 3 uniformly distributed beads.Considering more beads, uniformly distributed, does not affect the conformation tensor.
Fig. 4 compares two different configurations for r ¼ 5 and Fig. 5 depicts the angular velocity distribution for r ¼ 10.From these figures it can be noticed that despite the fact that the equivalent ellipsoid aspect ratio seems quite different, the angular velocity distributions are very close because k becomes very close to one when the aspect ratio is large enough (r > 5).In fact, for high aspect ratios, the angular velocity becomes quite insensitive to the aspect ratio itself.This means that inverse identification from finite element solutions could become inaccurate and then the noticed gap in predictions must be relativized.In any case, as soon as the actual aspect ratio becomes large enough (r > 5), the angular velocity distribution approaches the one characteristic of rods (r ! 1).

Conclusions and perspectives
In this paper, we have proposed an original multiscale description of dilute suspensions of rigid particles consisting of nine conceptual bricks: the first three focus on the particle itself (definition of conformation, evolution of conformation and induced stress); the next three on the population (description, evolution and contribution to the stress); and the last three on the microstructural macroscopic description, its evolution and the associated stress expression.
We have illustrated this micro-to-macro description in the case of dilute suspensions of rigid rods in a Newtonian fluid, where there exist numerous works focusing on each one on these bricks.This description being general, however, we propose to consider it for describing any kind of suspension of rigid particles, and in particular dilute suspensions of rigid clusters composed of rigid rods which we addressed here.
At the level of the particle, here the rigid cluster composed of rigid rods, we proved that the rotary velocity only depends on a symmetric tensor c with unit trace that can be considered as the appropriate conformation tensor for describing cluster kinematics.When the conformation tensor c has two identical eigenvalues, we proved that the resulting kinematics is that of an equivalent ellipsoid expressed by Jeffery's equation.To check the proposed model  it was applied numerically to different configurations of ellipsoidal clusters and then generalized to more complex shapes, for example cylindrical clusters.The predicted angular velocity was compared to the one obtained by existing models for cylindrical particles, where an excellent agreement was found for large and small cylinder aspect ratios (representing rods and disks respectively); small deviations were noticed when r % 1 that depend on the particular cluster configuration.Kinematics of rigid clusters with arbitrary shape were compared in the 2D case with direct numerical simulations in [1] and both predictions were in excellent agreement.
In the present paper, we have only addressed the first brick of the proposed approach to the multiscale modeling of suspensions involving rigid clusters composed of rigid rods.In way of perspectives for future work, we conclude with an overview of the complete nine-step approach to be followed: 1.The cluster conformation has been successfully described by tensor c given in Eq. ( 56).This choice was motivated by the fact that all clusters having the same conformation tensor c have the same kinematics as proved Section 3. 2. Cluster conformation evolution.Knowing the cluster kinematics x, we can obtain the expression of _ p in order to define the conformation evolution from _ c. 3. Cluster contribution to the stress.The force acting on each rod bead involved in the rigid cluster must be taken into account by invoking Kramers' rule.4. Population description.The population of rigid clusters at the mesoscopic scale can be represented by means of either a discrete or a continuous description.In the discrete framework, the population is described from different computational individuals, each one characterized by its conformation tensor c i ; i ¼ 1; . . .; N .Within the continuous framework, the suspension is characterized by the pdf W, which now depends on the physical coordinates (space x and time t) and the conformation coordinate c.Thus, the pdf reads Wðx; t; cÞ.It gives the fraction of clusters that at position x and time t have a conformation given by c. 5. Description of the population evolution.Within the discrete framework, the population evolution is obtained by integrating the evolution equation for each cluster conformation _ c i .Within the continuous framework, the pdf Wðx; t; cÞ evolves according to its associated Fokker-Planck equation, with the knowledge of _ c. 6. Contribution of cluster population to the stress requires adding the contribution of all computational clusters, within the discrete framework, or to integrate in conformation space when proceeding within the continuous framework.The remaining bricks (2)-( 9) constitute work in progress within our group.

Appendix A. On the cluster rotary velocity
The expression of the inverse in Eq. ( 55) deserves some additional comments.If for a while we define tensor m according to m ¼ I À c, as matrix m is 3 Â 3 it satisfies its characteristic polynomial (according to the Cayley-Hamilton theorem): By developing the expression of m À1 and taking into account that m ¼ I À c, we obtain ðI À cÞ À1 ¼ a 1 I þ a 2 c 2 ; ðA:4Þ where the coefficients a 1 and a 2 , that depend on the invariants of c, are given by a 1 ¼ Until now, we addressed Jeffery's solution for infinite aspect ratio ellipsoids, i.e., slender bodies or rods.Jeffery's solution for an ellipsoid immersed into a Newtonian fluid flow whose unperturbed velocity field is described by the velocity field v is given by [29]: where k depends on the ellipsoid aspect ratio r (the length to diameter ratio): and where for the sake of clarity we write the term D : ðp pÞ in the matrix algebra format, by using the equivalence D : ðp pÞ p T Á D Á p.We now prove that in order to represent finite aspect ratio 2D particles it suffices to consider a rigid system composed of two rods, aligned perpendicularly to each other, and having lengths 2L 1 and 2L 2 ; L 1 P L 2 , as depicted in Fig. 6, such that the parameter k in Jeffery's Eq. (B.1) results from r ðB:3Þ In the present configuration the forces applied on bead L 1 p 1 and L 2 p 2 are if a and b are first-order tensors, the single contraction Á reads ða Á bÞ ¼ a j b j (); if a and b are first-order tensors, the dyadic product reads ða bÞ jk ¼ a j b k ; if a and b are first-order tensors, the cross product Â reads ða Â bÞ j ¼ jmn a m b n , where jmn are the components of the Levi-Civita tensor (also known as permutation tensor); if a and b are respectively second and first-order tensors, the single contraction Á reads ða Á bÞ j ¼ a jm b m ; if a and b are second-order tensors, the single contraction Á reads ða Á bÞ jk ¼ a jm b mk ; if a and b are second-order tensors, the double contraction '':'' reads ða : bÞ ¼ a jk b kj .

ð47Þ
where wðp; LÞ is the angular distribution of beads located at distance L from the cluster center of gravity G, and CðLÞ is the fraction of beads located at that distance L. The normality condition reads Z S wðp; LÞ dp ¼ 1; 8L: ð48Þ By defining the conformation tensor c L related to the population of beads located at distance L with respect to the cluster center of gravity G as c L ¼ Z S p p wðp; LÞ dp; ð49Þ all sums in the previous expressions must be substituted by the corresponding integrals in the length and orientation domains, L and S respectively, weighted by the distribution function wðp; LÞ CðLÞ.

LL 2 L 2 c
By introducing the mean square length b and tensor c according to b ¼ Z CðLÞ dL; L CðLÞ dL;

Fig. 3
Fig. 3 illustrates the case of a simple shear flow v T ¼ ð0; _ cz; 0Þ

7 .Cc 8 .
The macroscopic description uses a coarser description based on the moments of the pdf.Here, the simplest choice consists of the first moment Cðx; tÞ defined as Cðx; tÞ ¼ Z Wðx; t; cÞ dc: ð59Þ Microstructural macroscopic evolution.In order to derive the time evolution of the first moment C, we should consider its time derivative and propose adequate closure relations.9.The moment-based stress requires consideration of the stress expression obtained in step (6) above.Here again, introduction of a suitable closure relation is required.

Table 1
[43]arison of the resulting ellipsoid aspect ratios with the ones predicted by Zhang et al.[43]for cylinders.