On the multi‐scale description of electrical conducting suspensions involving perfectly dispersed rods

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents 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. On the multi‐scale description of electrical conducting suspensions involving perfectly dispersed rods Marta Perez, Emmanuelle Abisset-Chavanne, Anais Barasinski, Francisco Chinesta, Amine Ammar, Roland Keunings

An important recent observation is that the flow-induced concentration of CNTs is not uniform [15][16][17][18][19]. Indeed, CNTs have the tendency to aggregate (in what follows we do not make distintion between agglomeration and aggregation, we use the last term because it was the one considered in our former works). As a result, all properties and mechanisms must be reformulated in the context of suspensions and networks involving clusters composed of CNTs, instead of considering a population of perfectly distributed isolated CNTs [5].
In our recent studies [20][21][22], we have proposed kinetic theory models to predict the kinematics and rheology of either rigid or deformable aggregates of CNTs in a Newtonian fluid matrix. In addition to providing a simple description of a rich microstructure, the proposed models are able to point out collective effects that are observed experimentally [23]. In these studies, we considered hydrodynamic effects only. Electrical mechanisms are discussed in the present paper.
In terms of modeling and simulation, two different approaches are usually adopted to take account of electrical effects. The simplest one considers a given network and proceeds to evaluate its direct current (DC) electrical properties. To that purpose, three main steps are followed: (1) generation of the composite's microstructure, (2) creation of an equivalent resistance network corresponding to this microstructure, and (3) calculation of this network in a continuous or discrete manner [24].
The above approach allows for a detailed analysis of the electrical properties for a given microstructure, which usually remains "frozen".
The second approach consists in predicting the network itself, and then carrying out the electrical analysis. In this case, it is usual to proceed at the mesoscopic scale using methods of Dissipative Particle Dynamics (DPD) for describing packed assemblies of oriented fibers suspended in a viscous medium [25]. Computer simulations are performed in order to explore how the aspect ratio and degree of fiber alignment affect the critical volume fraction percolation threshold required to achieve electrical conductivity. The fiber network impedance is assessed using Monte Carlo simulations after establishing the structural arrangement with DPD. Thus, these simulations allow one to predict the microstructure (CNT dispersion and aggregation), and, even though most such simulations do not consider the flow coupling, there are no major difficulties to include it as well. Computational micromechanics approaches for modeling effective conductivity in CNT-nanocomposites in general were addressed in [26,27].
The main limitations of these two common approaches are that (1) they concern a computational domain that is only representative of a small region of the whole process and part, and (2) they analyze a particular configuration, which implies many individual solutions in order to perform a valuable statistical treatment of the results.
In the present work, we propose an alternative approach to evaluating electrical properties in flowing suspensions of perfectly dispersed CNTs in a Newtonian fluid (the dispersion is ensured by assuming an appropriate functionalization). Starting from a microscopic description, we derive both mesoscopic and macroscopic descriptions. The main advantage of mesoscopic models is their ability to address systems of macroscopic size, while keeping track of the detailed physics through a number of conformational coordinates for describing the microstructure and its time evolution. At the mesoscopic scale, the microstructure is defined by means of the orientation distribution function that depends on physical space, time and CNT orientation. The moments of this distribution constitute a coarser description often used in macroscopic modeling, at the cost of compulsory closure approximations whose impact on the results is either ignored or unknown. Finally, the modeling of particle contacts is addressed as they determine the final functional properties, in particular electrical conduction. Different descriptors of rod contacts will be proposed and analyzed.

Orientation induced by the electric and flow fields
In this section, we first give the equation governing the orientation of a rod immersed in a Newtonian fluid of viscosity η in presence of an electric field ǫ(x, t) and a velocity field v(x, t). Then, the proposed model will be coarsened for describing a population of rods within the framework of kinetic theory. Finally, a macroscopic model will be derived.

Microscopic description
We consider a suspending medium consisting of a Newtonian fluid in which are suspended N rigid slender rods (e.g. CNTs) of length 2L. As a first approximation, the fiber presence and orientation are considered not to affect the flow kinematics defined by the velocity field v(x, t), with x ∈ Ω ∈ R 3 .
The microstructure is described at the microscopic scale by the unit vector defining the orientation of each rod, i.e. p i , i = 1, . . . , N. In absence of electric field, one fiber can be defined by p or −p, which implies a symmetry property for the orientation distribution function. When considering the electric field induced charges, however, that symmetry is broken and the orientation is defined univocally. We assume that p points from the negatively-charge bead to the positive one ( Fig. 1).
If the suspension is dilute enough, rod-rod interactions can be neglected and a micromechanical model can then be derived by considering a single generic rod whose orientation is defined by the unit vector p.
We thus consider the system illustrated in Fig. 1, that consists of a rod immersed into a fluid flow, with two beads at the extremities that have respectively an electrical charge +q and −q assumed for the sake of simplicity permanent and independent of the applied electric field. Hydrodynamic forces are also assumed to act on both beads. Thus, the resultant force acting on bead pL reads where ξ is the friction coefficient and E ≡ qǫ. It has been assumed that the velocity gradient is constant at the scale of the rod (first-gradient theory).
Obviously, if F is applied on bead pL, then for the opposite bead −pL the resultant force reads Neglecting inertia, the balance of linear and angular momenta yields [28] and This can be rewritten as: where ṗ E is the contribution of electrostatic forces to the rotary velocity and ṗ J represents the hydrodynamic contribution that coincides with the Jeffery solution for infinite aspect ratio ellipsoids [29].

Mesoscopic description
Because the rod population is very large, the description that we just proposed, despite of its conceptual simplicity, fails to address the situations usually encountered in practice. For this reason, coarser descriptions are preferred. The first plausible coarser description applies a zoom-out, wherein the rod individuality is lost in favour of a probability distribution function [30][31][32][33].

Fokker-Planck equation
In the case of rods, one can describe the microstructure at a certain point x and time t from the orientation distribution function Ψ (x, t, p) which gives the fraction of rods that at position x and time t are oriented in direction p. Obviously, the function Ψ satisfies the normality condition: where S is the surface of the unit ball that defines all possible rod orientations.
The balance equation that ensures conservation of probability reads: (1) F(pL) = E + ξ L(∇v · p −ṗ), S Ψ (x, t, p) dp = 1, ∀x, ∀t, For inertialess rods in first-gradient flows, we have ẋ = v(x, t) and the rod rotary velocity is given by Eq. (5): Equation (7), combined with Eq. (8), is known as a Fokker-Planck equation. It is a suitable compromise between the macroscopic scale that defines the overall process, and a finer microscopic description of the electric field induced orientation of each individual rod. The price to pay is the increase of the model dimensionality, since the orientation distribution is defined in a high-dimensional domain, i.e. (x, t, p) ∈ Ω × I × S.
In the case of dilute suspensions of CNTs, Brownian forces have a randomizing effect. These effects can be accurately described collectively with a diffusion term in the Fokker-Planck equation: where the diffusion coefficient D r quantifies the randomizing nature of rod-rod interactions and Brownian forces, that in semi-concentrated flowing systems was assumed scaling with the shear rate [34].

Parametric solutions of the Fokker-Planck equation
The Fokker-Planck equation (9) is defined in a multidimensional space involving the physical space x, the time t and the conformational coordinates associated to the rod orientation p, with ṗ given by Eq. (8).
We proposed a few years ago [35,36] a discretization technique based on the use of separated representations in order to ensure that the complexity scales linearly with the model dimensionality. This new approach is now known as Proper Generalized Decomposition (PGD). 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 Here, both the one-dimensional functions F j i and the appropriate number of terms N are unknown a priori. The interested reader can refer to the research papers [37][38][39][40][41][42][43] and the recent book [44].
One of the most appealing features of this technique is its ability of solving multidimensional models. Thus, with the PGD, the parameters of a physical model can be considered as extra-coordinates, such that by solving only once the resulting multidimensional model one has access to the general parametric solution that can be used to evaluate the impact of each parameter on the solution (that is, for performing sensitivity analyses) or to perform fast inverse identification [44][45][46].
In the case of the Fokker-Planck model (9)- (8), which only involves as coordinates the space x, the time t and the conformational coordinates p, its solution depends parametrically on the following parameters: (1) the diffusion coefficient D r , (2) the effective electric field E/ξ L and (3) the components of the velocity gradient ∇v Within the PGD framework, and assuming a uniform effective electric field Ẽ = E/ξ L, the parametric orientation distribution Ψ (x, t, p, D r ,Ẽ, ∇v) is written in a separated form for circumventing the curse of dimensionality associated with standard mesh-based discretization techniques. The simplest separated representation reads In 3D physical space, this formulation involves 18 dimensions. With the PGD, the computational complexity is associated with the solution of some 2D or 3D problems related to the calculation of the functions F , and a series of 1D problems for calculating the remaining functions involved in Eq. (12). The generic PGD solution procedure, that is, the constructor of the unknown low-dimensional functions involved in Eq. (12), is described in detail in the book [44].

Macroscopic description
Fokker-Planck based descriptions are rarely considered in industrial applications precisely because the computational complexity that high dimensionality induced by the use of conformation coordinates implies. For this reason, mesoscopic models are commonly coarsened one step further to obtain macroscopic models defined in standard physical domains, involving only space and time coordinates.
At the macroscopic scale, the orientation distribution function is substituted by its moments for describing the microstructure, as proposed in [47]. We consider the first four orientation moments a (1) , a (2) , a (3) and a (4) defined as: In view of the non-symmetry of Ψ, the odd moments of the distribution function Ψ do not vanish.

Evolution equations for the orientation moments
As detailed in [28], by taking the time derivative of Eqs. (13) and (14) and taking into account the expression of the rotary velocity ṗ expressed by Eq. (5), we obtain and both requiring closure relations, the first one (17) concerning a (2) and a (3) and the second one (18) associated with a (3) and a (4) .
When Brownian effects are retained, the resulting macroscopic model is given by [28] which, in absence of electric field, E = 0 and of flow, v(x, t) = 0, ensures a steady-state isotropic distribution, i.e. a (1) The second moment evolves according to where d = 2 in 2D and d = 3 in 3D. In absence of electric field and of flow, we obtain a steady-state isotropic distribution, i.e. a (2)

On closure relations
We proposed and analyzed in [28] the following hybrid closure relation for a (2) , to be considered when the microstructure evolution is described from the single a (1) descriptor requiring the integration of Eq. (17): with d = 2, 3 in 2D and 3D respectively, and Integration of Eq. (17) also requires the third-order moment. In [28], the following pragmatic choice was considered even though, as discussed below, it does not satisfy the required symmetry conditions. When considering both microstructure descriptors, a (1) and a (2) , closure relations are needed for approximating the third and fourth-order moments, a (3) and a (4) respectively.
There is a vast choice of possible closure relations of the fourth-order tensor a (4) usually encountered in standard suspension models [48][49][50].
When there is no difference between p and −p to describe rod orientation, the distribution function is symmetric and odd moments vanish. In the present case, the symmetry is broken and odd moments do not vanish, thus requiring appropriate closures.
From the definition of a (3) , we can identify the following symmetry conditions (for the sake on simplicity and without loss of generality we restrict the analysis to the 2D case): We could consider the following terms in the closure of the third-order moment: -The cubic term a (1) ⊗ a (1) ⊗ a (1) fulfills the symmetry relations (25), whereas a (1) ⊗ a (2) and a (2) ⊗ a (1) do not satisfy in general conditions (25); -Quadratic terms obtained from tensor products of a (2) and constant vectors or quadratic terms coming from the tensor product of a (1) twice and constant vectors do not verify in the general case the symmetry conditions; -Linear terms obtained from tensor products of a (1) and constant matrix do not verify in the general case the symmetry conditions; -If we define vectors I T 1 = (1, 0) and I T 2 = (0, 1), then tensors I 1 ⊗ I 1 ⊗ I 1 and I 2 ⊗ I 2 ⊗ I 2 verify the above mentioned symmetry conditions.
Thus, we could consider the following closure: In the numerical analysis and discussion reported later, we will consider the cubic closure (26) with β = δ = 0, i.e.

Description of the rod network
We have seen how to model and predict the microstructure induced by the electric field. We addressed in [28] the following questions: (1) how to compute the electric field (1) .
E(x, t), (2) how to evaluate the induced conductivity properties, and finally (3) how to determine preferential electrical paths in the computational domain of interest Ω.
In what follows, the electric field is assumed known from the solution of the Laplace equation in the domain occupied by the suspension (see [28] for details). We focus here on the multi-scale description of rod contacts.

Mesoscopic description
In order to quantify the rod network, we introduce the number density of rod contacts C(x, t, p) for a rod with orientation p, depending on the two main microstructure descriptors: (1) the CNT concentration φ(x, t), and (2) the orientation distribution Ψ (x, t, p).
When calculating the number of contacts, there is no difference between p and −p. Thus, we use the symmetrized orientation distribution function Ψ S (x, t, p) defined as follows: The fraction of rods oriented in direction p is φ Ψ S (p). The center of gravity of all of these is located at position P. All rods having their center of gravity at position Q inside the sphere of radius 2L centered at point P could interact with the former. Due to the small size of the considered rods (CNTs), we can assume that the orientation distribution Ψ and concentration φ at positions P and Q coincide.
In fact, because of electronic tunneling effect, the contact between the fibers is not necessarily needed. There is a minimum distance, δ, from which the electrical conductivity occurs. The interaction between rods aligned in directions p and p ′ whose centers of gravity are closer than 2L can be evaluated by a simple geometrical construction as shown in Fig. 2.
The local and directional contact density C(p), when its dependence on the space and time coordinates is omitted, reads (29) C(p) = φ B S χ(r, p, p ′ ) Ψ S (p ′ ) dp ′ dr,

Fig. 2 Evaluating rod contacts
where B is the ball of radius 2L centered at point P and where χ(r, p, p ′ ) is a function equal either to 1 if the distance between the fibers is lower than δ, or to zero otherwise.
It can be expected that the conductivity at position x and time t along the direction p scales with the number of contacts φ C(x, t, p)Ψ S (x, t, p). Thus the conductivity becomes local and directional. One can expect that percolation along direction p occurs locally (at position x and time t) when the number of contacts φ C(x, t, p)Ψ S (x, t, p) is higher than a threshold value T . Lower values imply an infinite local directional electrical resistivity. Thus, percolation could be considered local and directional, allowing us to create a network as described in [28] in which we associate to each node a connexion with each one of its 26 neighbours (Fig. 3). Once we know the local directional resistivity at node (i, j, k) scaling with φ C(x, t, p)Ψ S (x, t, p), we can compute the electrical resistances between node (i, j, k) and each one of its neighbouring nodes, and, using Kirchhoff 's law, analyze the resulting electrical circuit as described in [28].
The microstructure description based on the use of the density of contacts C(p) is very close to that obtained following the rationale proposed by Toll [51,52] (in continuity with the works of Doi and Edwards [53] and Ranganathan and Advani [54]). In our approach, the effect related to the finite diameter of the rods is neglected, but its inclusion is straightforward.
The integration involved in Eq. (29) can be understood from the geometrical construction depicted in Fig. 4. It can be noticed that all rods oriented along direction p ′ having their centre of gravity inside the parallelogram of side length 2L, defined by the black broken line, intersect the test rod (the horizontal one). The area of this parallelogram is (2L) 2 �p × p ′ �. Thus, integral (29) can be rewritten as The density of contacts C(p) is appropriate because it constitutes a rich description of the microstructure. Its calculation, however, requires knowledge of the distribution function Ψ (p) that is only available when operating at the mesoscale by solving the Fokker-Planck equation.

Towards a fully-macroscopic description
The question arising immediately concerns the existence of an appropriate fully-macroscopic description.
The simplest proposal consists in defining a tensor L as follows:

which in view of Eq. (30) can be rewritten as
This expression is very close to the interaction tensor b introduced by Ferec et al. [55], Indeed, we have the equality As the distribution function Ψ S (p) can be approximated from its associated moments (only the even ones in view of symmetry), as proven in [47], the interaction tensor b can be written from the different orientation tensors a S,(n) , n = 2, 4, . . ., where the superscript S highlights the fact that all of them are associated to the symmetric orientation distribution function Ψ S (p). It is easy to verify from the moments' definition that This procedure for expressing b as a function of a (2) and a (4) was considered in [56] where an appropriate expression b = b(a (2) , a (4) ) was proposed and validated. This route allows us to finally establish a link between our descriptor L (32) and the standard orientation tensors, which allows for a fully-macroscopic description. Now, from the knowledge of L, we can obtain the directional properties along direction p by calculating (p T · L · p) p, or equivalently (L : (p ⊗ p)) p.
Another appealing microstructural descriptor could be the second-order tensor J defined such as to verify In general, such a tensor does not exist, and then the most natural possibility is to enforce this equality is a least-squares sense. Thus, we define the functional F(J),   This expression constitutes a sort of bridge (in a least-squares sense) between macroscopic and mesoscopic scales. Indeed, knowing L (which can be evaluated from the orientation tensors a (2) and a (4) ) and using Eq. (39), we can evaluate the tensor J, from which we calculate C(p) (in a least-squares sense) from (36).
It is important to recall that, in the least-squares approximation (37), the equality is weighted by the distribution function Ψ S (p). This means that more than reproducing C(p), we are approximating C(p)Ψ S (p), which constitutes a better description from the physical point of view, as it combines the potential of a given direction expressed by C(p) with the number Ψ (p) of rods having this orientation.

Evaluating the closure relations
We now evaluate and validate the different closure relations introduced above.
In "Macroscopic description", we derived two frameworks for describing the microstructure evolution. The first one considers the first moment a (1) such that It requires two closure relations for expressing the higher-order moments a (2) and a (3) as a function of a (1) .
The second framework involves the first two moments, where Eq. (40) is complemented with an equation for the second-order moment: This second framework requires appropriate closures for the third and fourth-order moments.
First, for different orientation distributions Ψ (p, t) coming from the solution of the Fokker-Planck equation for different choices of the velocity and electric fields, we compute the third-order moment a (3),FP from Eq. (15) (the superscript FP refers to the Fokker-Planck nature of the solution). We then fit the coefficient α in Eq. (27) in order to minimize the gap (in a least-squares sense) between a (3),FP and a (3),cl , The best fit was obtained with α adjusted from several numerical experiments as follows: It is easy to verify by taking the trace of both sides of Eq. (41) that d dt tr a (2) (t) = 0, which implies tr(a (2) (t)) = 1 if the initial condition has a unit trace. However, as the closure relation (45,46) does not reproduce exactly a (3) , we must check the model consistency by examining if the trace of a (2) (t) remains unitary when integrating Eq. (41) from an initial condition with unit trace.
In simple shear flow, integration of Eq. (41) with the closure relation (45,46) yields an unphysical variation of the trace of a (2) . This fact can be observed in Fig. 5 which depicts the evolution of a (2) (t) (on the left) and its trace (on the right) resulting from the integration of Eq. (41) for planar shear flow v T = (y, 0), Ẽ = 1, D r = 0.1 and an isotropic planar orientation state.
In relation to Eq. (41), the condition that the third-order closure must satisfy in order to keep the trace constant is 2a (1) · E = 2 tr a (3) · E that can be easily verified (43) a (4),qua = a (2) ⊗ a (2) ij I kl + a (2) ik I jl + a (2) il I jk + a (2) kl I ij + a (2) jl I ik + a (2) jk I il , (1) . When considering the closure (45) into Eq. (41), we obtain that must be equal to 2a (1) · E in order to guarantee a constant value of the trace of the second-order moment. The only possibility is that α in Eq. (45) verifies from which the third-order moment closure verifying both the symmetry conditions and the model consistency (with respect to the preservation of the unit trace of a (2) ) must be The closure relation (49), however, when considered in the integration of Eq. (41), yields large deviations with respect to the reference solution obtained from the Fokker-Planck solution, as can be seen in Fig. 6 for the same flow conditions as considered previously.

Fig. 6
Second order moment a (2) (t) computed with the closure (50) from which the trace of Eq. (41) involving the electric contribution gives where E 1 and E 2 are the components of E. Considering β = νa (1) 1 and δ = υa (1) 2 , the previous expression can be rewritten as which must be equal to 2a (1) · E in order to guarantee a constant value of the trace of the second-order moment. This equality holds if µ + ν + υ = 1. All the possibilities that such richer closures open constitute work in progress whose results will be reported in further works. In what follows, we consider an alternative route.
This alternative route consists in reinterpreting the pragmatic closure (24) by assuming that it is in fact a closure of the tensor product a (3) · E, that is With the second moment closure verifying the normality property tr a (2),cl = 1, as ensured for example by Eq. (21), we have the equality that ensures the model consistency (unit trace).
Thus, in what follows, the retained closures of the third-order terms in Eqs. (40) and (41) are respectively and In the forthcoming numerical examples, we use Eqs. (40) and (41) with the fourth-order hybrid closure (42) and the third-order closures (57) and (58). We first consider planar flows characterized by a velocity gradient of the general form and different values of the diffusion coefficient. (1) tr(a (1) ⊗ a (1) ) + β(I 1 ⊗ I 1 ⊗ I 1 ) + δ(I 2 ⊗ I 2 ⊗ I 2 ), In all cases, the Fokker-Planck equation (9) is solved for the orientation distribution Ψ (p, t) by means of the PGD. The initial orientation distribution is assumed isotropic. From Ψ (p, t), the first and second-order moments a (1),FP (t) and a (2),FP (t) are calculated, respectively from Eqs. (13) and (14). This procedure does not involve any closure relation. The solutions a (1),FP (t) and a (2),FP (t) can thus be considered as reference solutions.
The same moments are now computed via Eqs. (40) and (41) using closure approximations. When integrating Eq. (40) to obtain a (1) (t), we take the second-order moment from the solution of Eq. (41) and the third-order moment is approximated by using the closure (57). When integrating Eq. (41) to obtain a (2) (t), we use the closure (58) for the third-order moment and the standard hybrid closure (42) for the fourth-order moment.
We consider both simple shear and extensional flows for evaluating the accuracy of the resulting closed macroscopic description: In Figs. 7 and 8, the approximate closure solutions a (1) and a (2) are compared with the reference solutions a (1),FP and a (2),FP obtained from the direct solution of the Fokker-Planck equation. Agreement is in general quite satisfactory, but the development of even more precise closures will require a considerable amount of work.
It is well known that closure relations can induce artifacts, as for example the disappearance of hysteric behaviour in reversal flows. To check the model response in such conditions, we apply a simple shear flow defined by v T = (γ y, 0, 0) for t ∈ [0, T = 5s]; the flow is reversed for t ∈ [T , 2T ] by changing the sign of γ, and it is again reversed in subsequent time intervals. In other words, the shear flow is reversed at times t = nT, n = 1, 2, . . . We see in Fig. 9 that the closure solution a (1) constitutes a reasonable approximation of the reference solution, again computed from the orientation distribution solution of the associated Fokker-Planck equation. In order to appreciate this more clearly, Fig. 10 compares the component a (1) y of both solutions. Finally, we compare in Fig. 11 the planar component of a (2) and a (2),FP . All these results show that the proposed closure relations remain quite reasonable and do not introduce unphysical behaviours.

Parametric solution of the Fokker-Planck equation
When controlling the applied flow and electric fields, one could require many solutions of the model in order to reach optimal conditions with respect to an output of interest. In these circumstances, the parametric solution of the model would be extremely In what follows, we consider a homogeneous simple shear flow v T = (γ y, 0, 0), γ = 1. Its homogeneity avoids the dependence of the orientation distribution function Ψ on the space coordinates x. We assume planar orientation, with a single conformational coordinate to describe it, i.e. the angle θ with respect the x−coordinate axis. An electric field of intensity E is applied along the y−coordinate axis, i.e. E = Ei y (with i y the unit vector defining the y−coordinate axis direction). In this section, E refers to the effective electric field denoted by Ẽ in "Parametric solutions of the Fokker-Planck equation".
Thus, the four-dimensional parametric solution Ψ (t, θ , D r , E) of the Fokker-Planck equation (7) is sought in the separated form  Figure 12 depicts the different functions involved in the separated representation (60). Finally, Fig. 13 shows the reconstructed parametric steady state solution Ψ (t → ∞, θ , D r , E). For more details on the practical implementation of the PGD for solving parametric PDE's, see [44][45][46].

Evaluating interaction descriptors
Finally, we check whether the different interaction descriptors discussed in "Description of the rod network" perform as expected.
First, given the two planar orientation distributions shown in Fig. 14, we compute the density of contacts C(θ) given by Eq. (29) and the one given by Toll [51,52]. Both results are compared in Fig. 15. We observe an excellent agreement despite of some hypotheses introduced in our approach (related to the assumption of a vanishing rod diameter). As expected, Fig. 16 shows the perfect agreement between the scaled Ferec interaction tensor b, Eq. (33), and the one that we proposed above, L, Eq. (31).
Finally, we discussed in "Towards a fully-macroscopic description" the use of tensor J as an interaction descriptor in order to fit in a least-squares sense the density of contacts. The least-squares being weighted by the distribution function Ψ, the resulting directional information J : (p ⊗ p) is closer to CΨ than to C, as can be noticed in Fig. 17.

Conclusions
In this work, we have revisited the microstructural kinematics of electrically conductive rods immersed in a flowing suspension to which an electric field is applied. The evolution equation for the microscopic rod orientation was derived by considering standard balances of forces and moments. These kinematics were then introduced into a mesoscopic description based on the Fokker-Planck equation, whose parametric solution was formulated within the separated representation framework at the heart of the method of Proper Generalized Decomposition. The mesoscopic description was coarsened one step further in order to derive a fully macroscopic description based on the use of the first or the first two moments of the orientation distribution. Both approaches require the use of appropriate closure relations. In the present work, we proposed consistent closures for the third-order orientation tensor, which despite being non optimal produce results that are in reasonable agreement with reference solutions obtained directly from the Fokker-Planck description.
On the other hand, prediction of the electrical properties of the suspension as a whole requires the quantification of the rod interaction. We quantified the directional conductivity by calculating the directional density of contacts that depends on the orientation distribution and the local rod concentration. This approach, already considered in our former work [28], requires however the mesoscopic calculation of the orientation distribution. In this work, we proposed an alternative fully-macroscopic route. For that purpose, we defined an interaction tensor L which by construction is equivalent to the interaction tensor proposed by Ferec in [55]. As the latter can be approximated with reasonable accuracy from the orientation tensors of even order, a fully macroscopic approach is thus finally obtained.