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

Emmanuelle Abisset-Chavanne, Francisco Chinesta, Julien Férec, Gilles Ausias, Roland Keunings

► To cite this version:

Emmanuelle Abisset-Chavanne, Francisco Chinesta, Julien Férec, Gilles Ausias, Roland Keunings. On the multiscale description of dilute suspensions of non-Brownian rigid clusters composed of rods. Journal of Non-Newtonian Fluid Mechanics, 2015, 222, pp.34-44. 10.1016/j.jnnfm.2014.08.014. hal-01584797

HAL Id: hal-01584797

https://hal.science/hal-01584797

Submitted on 9 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
On the multiscale description of dilute suspensions of non-Brownian rigid clusters composed of rods

E. Abisset-Chavanne<sup>a</sup>, F. Chinesta<sup>a</sup>, J. Férec<sup>b</sup>, G. Ausias<sup>b</sup>, R. Keunings<sup>c</sup>

<sup>a</sup> GeM, UMR CNRS-Centrale Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France

<sup>b</sup> LIMATB, Université de Bretagne-Sud, rue de St Maudé, F-56325 Lorient, France

<sup>c</sup> ICTEAM, Université catholique de Louvain, Bat. Euler, Av. Georges Lemaitre 4, B-1348 Louvain-la-Neuve, Belgium

The motion of an ellipsoidal particle immersed in a flow of a Newtonian fluid was obtained in the pioneering 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 flow, 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 flow of a Newtonian fluid. 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.

Keywords:

Aggregate

Cluster

Jeffery's equation

Non-ellipsoidal particles

Multiscale modeling

Mesoscopic kinetic theory

# 1. 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 –pdf– whose 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 be 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, \dots, x_d)$ as:

$$u(x_1, \dots, x_d) \approx \sum_{i=1}^{i=N} F_i^1(x_1) \dots F_i^d(x_d). \qquad (1)$$

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–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 $\dot{\mathbf{p}}$ of the particle conformation. For an ellipsoid suspended in a Newtonian viscous fluid, whose kinematics is described by the gradient of velocity $\nabla \mathbf{v}$, one obtains Jeffery's equation [29]:

$$\dot{\mathbf{p}} = \mathbf{\Omega} \cdot \mathbf{p} + \lambda(\mathbf{D} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}), \qquad (2)$$

where $\mathbf{\Omega} = \frac{1}{2}(\nabla \mathbf{v} - (\nabla \mathbf{v})^T)$ and $\mathbf{D} = \frac{1}{2}(\nabla \mathbf{v} + (\nabla \mathbf{v})^T)$ are respectively the skew-symmetric and symmetric components of the gradient of velocity tensor, and $\lambda$ is a coefficient that depends on the ellipsoid aspect ratio $r$ (ratio of length $2L$ to diameter $\phi$): $\lambda = \frac{r^2-1}{r^2+1}$. In view of the symmetry of the strain rate tensor $\mathbf{D}$, the equality $\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}) = \mathbf{D} : (\mathbf{p} \otimes \mathbf{p})$ applies.

Three special cases can be considered:

* Slender body theories result for infinite aspect ratio ellipsoids, implying $r \to \infty$ and then $\lambda \to 1$. Thus Jeffery's equation reduces to:

$$\dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} + (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}. \qquad (3)$$

* When the ellipsoid length $2L$ coincides with the diameter of its circular cross section $\phi$ the ellipsoid reduces to a sphere and in that case $r = 1$ and $\lambda = 0$, implying that the sphere rotates according to the fluid flow vorticity $\mathbf{\Omega}$, that is:

$$\dot{\mathbf{p}} = \mathbf{\Omega} \cdot \mathbf{p}. \qquad (4)$$

* 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 \to 0$ and $\lambda \to -1$. The disk kinematics is then given by:

$$\dot{\mathbf{p}} = \mathbf{\Omega} \cdot \mathbf{p} - \mathbf{D} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}, \qquad (5)$$

where $\mathbf{p}$ is the unit vector normal to the disk surface.

5. **Particle contribution to the stress.**

4. **Population description.** Description of a population of $N$ ellipsoidal particles located at position $\mathbf{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 $\mathbf{p}_i$, $i = 1, \dots, N$ (discrete description); or from the introduction of a probability distribution function $-\text{pdf-} \psi(\mathbf{x}, t, \mathbf{p})$ that gives the fraction of particles that at position $\mathbf{x}$ and time $t$ are aligned along the direction $\mathbf{p}$ (continuous description). In the case of homogeneous distributions, the space dependence can be ignored and then $\psi(t, \mathbf{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 $\dot{\mathbf{p}}_i$ according to Jeffery's equation for each particle $i = 1, \dots, N$, or from the integration of the Fokker-Planck equation governing the evolution of the pdf $\psi(\mathbf{x}, t, \mathbf{p})$. This equation expresses conservation of probability. It reads:

$$\frac{\partial \psi}{\partial t} + \nabla_x \cdot (\dot{\mathbf{x}} \psi) + \nabla_p \cdot (\dot{\mathbf{p}} \psi) = 0, \qquad (6)$$

where $\nabla_x$ and $\nabla_p$ refer respectively to the gradients in spatial and conformational coordinates. If the vector $\mathbf{p}$ is defined from the ellipsoid center of gravity $\mathcal{G}$ and inertia effects are neglected, $\dot{\mathbf{x}} = \mathbf{v}$ as it can be shown that the particle center of gravity $\mathcal{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 $(\psi(\mathbf{p}) = \psi(-\mathbf{p}))$ in the case of rigid rods addressed here, odd-order moments vanish, and the first non-vanishing moment is the so-called second-order orientation tensor defined as:

$$\mathbf{a} = \int_S \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}) \, d\mathbf{p}, \qquad (7)$$

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:

$$ \mathbf{A} = \int_S \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}) \, d\mathbf{p}. \qquad (8) $$

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 sixth-order 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:

* if **a** and **b** are first-order tensors, the single contraction · reads ($\mathbf{a} \cdot \mathbf{b}$) = $a_j b_j$ ();

* if **a** and **b** are first-order tensors, the dyadic product ⊗ reads ($\mathbf{a} \otimes \mathbf{b}$)$_{jk}$ = $a_j b_k$;

* if **a** and **b** are first-order tensors, the cross product × reads ($\mathbf{a} \times \mathbf{b}$)$_j$ = $\epsilon_{jmn} a_m b_n$, where $\epsilon_{jmn}$ are the components of the Levi-Civita tensor $\epsilon$ (also known as permutation tensor);

* if **a** and **b** are respectively second and first-order tensors, the single contraction · reads ($\mathbf{a} \cdot \mathbf{b}$)$_j$ = $a_{jm} b_m$;

* if **a** and **b** are second-order tensors, the single contraction · reads ($\mathbf{a} \cdot \mathbf{b}$)$_{jk}$ = $a_{jm} b_{mk}$;

* if **a** and **b** are second-order tensors, the double contraction ":" reads ($\mathbf{a} : \mathbf{b}$) = $a_{jk} b_{kj}$.

## 2. Multiscale modeling of dilute suspensions of non-Brownian rods

We consider a suspending medium consisting of a Newtonian fluid of viscosity $\eta$ 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 $\mathbf{v}(\mathbf{x}, t)$, with $\mathbf{x} \in \Omega \in \mathbb{R}^d$, $d = 2, 3$. In what follows we consider in detail the different conceptual bricks previously introduced.

### 2.1. 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 $\mathcal{G}$ and aligned along the rods axis. Inertial effects are neglected in the sequel.

### 2.2. 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 $\mathbf{v}_0 + \nabla\mathbf{v} \cdot \mathbf{p}L$ and the second one by $\mathbf{v}_G + \dot{\mathbf{p}}L$. Thus, the force $\mathbf{F}(\mathbf{p}L)$ reads:

$$ \mathbf{F}(\mathbf{p}L) = \xi(\mathbf{v}_0 + \nabla\mathbf{v} \cdot \mathbf{p}L - \mathbf{v}_G - \dot{\mathbf{p}}L), \qquad (9) $$

where $\xi$ is the friction coefficient, $\mathbf{v}_0$ the fluid velocity at the rod center of gravity, and $\mathbf{v}_G$ the velocity of the center of gravity.

If **F** applies on the bead at $\mathbf{p}L$, then the force on the opposite bead at $-\mathbf{p}L$ reads

$$ \mathbf{F}(-\mathbf{p}L) = \xi(\mathbf{v}_0 - \nabla\mathbf{v} \cdot \mathbf{p}L - \mathbf{v}_G + \dot{\mathbf{p}}L). \qquad (10) $$

By adding Eqs. (9) and (10), and neglecting inertial effects, we obtain the force balance

$$ \mathbf{F}(\mathbf{p}L) + \mathbf{F}(-\mathbf{p}L) = 2\xi(\mathbf{v}_0 - \mathbf{v}_G) = \mathbf{0}, \qquad (11) $$

which implies $\mathbf{v}_0 = \mathbf{v}_G$, that is, the rod center of gravity is moving with the fluid velocity. For simplicity of notation, we shall write $\mathbf{F} = \mathbf{F}(\mathbf{p}L)$ and $\mathbf{F}(-\mathbf{p}L) = -\mathbf{F}$.

As the resulting torque must also vanish, the only possibility is that force **F** acts along **p**, that is $\mathbf{F} = \lambda\mathbf{p}$, with $\lambda \in \mathbb{R}$. Thus, we can write

$$ \lambda\mathbf{p} = \xi L(\nabla\mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}). \qquad (12) $$

Premultiplying Eq. (12) by **p** and taking into account that $\mathbf{p} \cdot \mathbf{p} = 1$ and consequently $\mathbf{p} \cdot \dot{\mathbf{p}} = 0$, we have

$$ \lambda = \xi L(\nabla\mathbf{v} : (\mathbf{p} \otimes \mathbf{p})), \qquad (13) $$

which gives

$$ \xi L(\nabla\mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p} = \xi L(\nabla\mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}). \qquad (14) $$

We have thus obtained Jeffery's Eq. (3),

$$ \dot{\mathbf{p}} = \nabla\mathbf{v} \cdot \mathbf{p} - (\nabla\mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p}. \qquad (15) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Hydrodynamic forces applied on a rod immersed in a Newtonian fluid.</figcaption></figure>
**Remark 2.** As the factor $\xi L$ appears in both sides of Eq. (14), the rod kinematics does not contain size effects.

### 2.3. Particle contribution to the stress

The forces applied at the rod ends $\mathbf{p}L$ and $-\mathbf{p}L$ are respectively $\lambda \mathbf{p}$ and $-\lambda \mathbf{p}$, i.e. directed along the rod and in equilibrium by construction.

With $\lambda$ given by Eq. (13), we have

$$ \mathbf{F}(\mathbf{p}L) = \xi L(\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}. \quad (16) $$

By applying Kramers' formula, the corresponding contribution to the stress is given by

$$ \boldsymbol{\tau}^p = \xi L^2 (\nabla \boldsymbol{v} : (\boldsymbol{p} \otimes \boldsymbol{p})) \boldsymbol{p} \otimes \boldsymbol{p}, \quad (17) $$

which can be rewritten as

$$ \boldsymbol{\tau}^p = \xi L^2 \nabla \boldsymbol{v} : (\boldsymbol{p} \otimes \boldsymbol{p} \otimes \boldsymbol{p} \otimes \boldsymbol{p}). \quad (18) $$

**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 $\mathbf{T}$ applied on the rod end. Considering Newtonian behavior, we have $\mathbf{T} = 2\eta \mathbf{D} \cdot \mathbf{p}$. By writing the objective stress related to the rod rotation $2\eta(\dot{\mathbf{p}} - \mathbf{\Omega} \cdot \mathbf{p})$ and enforcing that the resulting force aligns along the direction $\mathbf{p}$, we obtain again Jeffery's equation. In order to obtain the same expression for the stress, it suffices to consider $\xi \propto \eta/L^2$. Thus, we can conclude that size effects are absent in Eq. (18).

### 2.4. 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 $\mathbf{p}_i$, $i = 1, \dots, 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 $\psi(\mathbf{x}, t, \mathbf{p})$ that gives the fraction of rods that a position $\mathbf{x}$ and time $t$ are oriented along direction $\mathbf{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. $\mathbf{p}$ in the present case.

### 2.5. Description of the population evolution

* When the population is described from the individuals composing it, whose conformation is given by vectors $\mathbf{p}_i$, $i = 1, \dots, N$, the evolution of each one is given by Jeffery's equation:

$$ \dot{\mathbf{p}}_i = \nabla \mathbf{v} \cdot \mathbf{p}_i + (\nabla \mathbf{v} : (\mathbf{p}_i \otimes \mathbf{p}_i)) \mathbf{p}_i, \quad \forall i = 1, \dots, N. \quad (19) $$

* The alternative description consists of using the pdf $\psi(\mathbf{x}, t, \mathbf{p})$ that satisfies the normalization condition:

$$ \int_S \psi(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p} = 1, \quad \forall \mathbf{x}, \quad \forall t. \quad (20) $$

Conservation of probability yields

$$ \frac{\partial \psi}{\partial t} + \nabla_{\mathbf{x}} \cdot (\dot{\mathbf{x}} \psi) + \nabla_{\mathbf{p}} \cdot (\dot{\mathbf{p}} \psi) = 0, \quad (21) $$

where, for inertialess rods, $\dot{\mathbf{x}} = \mathbf{v}(\mathbf{x}, t)$, and the rod rotary velocity is given by Jeffery's equation:

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p}. \quad (22) $$

The price to pay is the increase of the model dimensionality, as the orientation distribution is defined in a 6-dimensional domain, i.e. $\psi : (\mathbf{x}, t, \mathbf{p}) \to \mathbb{R}^+$ where $\mathbf{x} \in \Omega \subset \mathbb{R}^d, d = 3; t \in \mathcal{I} \subset \mathbb{R}^+, \mathbf{p} \in S$.

### 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 $\mathbf{p}_i$ vectors, the contribution of rods to the suspension stress is calculated by adding their individual effects given by Eq. (18), that is:

$$ \boldsymbol{\tau}(\mathbf{x}, t) = \sum_{i=1}^{N(\mathbf{x}, t)} \boldsymbol{\tau}^i = \sum_{i=1}^{N(\mathbf{x}, t)} \xi L^2 \nabla \boldsymbol{v} : (\mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i), \quad (23) $$

where $N(\mathbf{x}, t)$ refers to the number of computational rods located at time $t$ in the neighborhood of $\mathbf{x}$.

* When the population is described from the pdf, the sum in Eq. (23) is replaced by an integral in the conformation space $S$:

$$ \begin{aligned} \boldsymbol{\tau}(\mathbf{x}, t) &= \int_S \boldsymbol{\tau}^p \psi(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p} \\ &= 2\eta N_p \int_S \nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}) \psi(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p}. \end{aligned} \quad (24) $$

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:

$$ \boldsymbol{\tau} = 2\eta N_p (\mathbf{A} : \nabla \mathbf{v}), \quad (25) $$

which in view of the symmetry of $\mathbf{A}$ can be rewritten as

$$ \boldsymbol{\tau} = 2\eta N_p (\mathbf{A} : \mathbf{D}). \quad (26) $$

### 2.7. 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.

## 2.8. 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:

$$ 
\begin{aligned}
\dot{\mathbf{a}} &= \int_S (\dot{\mathbf{p}} \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}) \psi \, d\mathbf{p} \\
&= \int_S (\nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}) \otimes \mathbf{p} \psi \, d\mathbf{p} \\
&\quad + \int_S \mathbf{p} \otimes (\nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}) \psi \, d\mathbf{p} \\
&= \nabla \mathbf{v} \cdot \mathbf{a} + \mathbf{a} \cdot (\nabla \mathbf{v})^T - 2 \mathbf{A} : \nabla \mathbf{v}.
\end{aligned}
\quad (27) $$

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 lower-order 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:

$$ \mathbf{A} \approx \mathbf{a} \otimes \mathbf{a}. \quad (28) $$

This gives

$$ \dot{\mathbf{a}} \approx \nabla \mathbf{v} \cdot \mathbf{a} + \mathbf{a} \cdot (\nabla \mathbf{v})^T - 2 (\nabla \mathbf{v} : \mathbf{a}) \mathbf{a}, \quad (29) $$

and invoking again symmetry considerations,

$$ \dot{\mathbf{a}} \approx \nabla \mathbf{v} \cdot \mathbf{a} + \mathbf{a} \cdot (\nabla \mathbf{v})^T - 2 (\mathbf{D} : \mathbf{a}) \mathbf{a}. \quad (30) $$

## 2.9. Moment based stress

We obtained previously the expression of the rod population contribution to the stress:

$$ \tau = 2\eta N_p (\mathbf{A} : \nabla \mathbf{v}), \quad (31) $$

which involves the fourth-order moment A. There is no closure issues when A is calculated from the pdf $\psi$ 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

$$ \tau = 2\eta N_p (\mathbf{A}^{cr}(\mathbf{a}) : \nabla \mathbf{v}), \quad (32) $$

where the superscript cr refers to the use of an appropriate closure relation.

With the quadratic closure, the stress reads:

$$ \tau = 2\eta N_p (\mathbf{a} : \nabla \mathbf{v}) \mathbf{a} \equiv 2\eta N_p (\mathbf{a} : \mathbf{D}) \mathbf{a}. \quad (33) $$

# 3. 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–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.

## 3.1. Discrete description

First, we consider a 3D rigid cluster consisting of N/2 rods $\mathcal{R}_j$ of length $Q_j$. We assume that each rod $\mathcal{R}_j$ contains two beads at its ends on which hydrodynamic forces apply. Thus, the cluster contains N beads $B_i$, $i = 1, \dots, N$. The location of each bead $B_i$ with respect to the cluster center of gravity $\mathcal{G}$ is represented by $L_i \mathbf{p}_i$, where $\mathbf{p}_i$ is the unit vector pointing from $\mathcal{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 $\mathbf{F}_i$ apply on each bead located at position $L_i \mathbf{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:

$$ \mathbf{F}_i = \xi (\mathbf{v}_0 + \nabla \mathbf{v} \cdot \mathbf{p}_i L_i - \mathbf{v}_G - \dot{\mathbf{p}}_i L_i), \quad (34) $$

where $\xi$ is the friction coefficient, $\mathbf{v}$ the flow velocity field, $\mathbf{v}_0$ the fluid velocity at position $\mathcal{G}$ and $\mathbf{v}_G$ the velocity of the center of gravity.

By adding all the forces we obtain

$$ \mathbf{0} = \sum_{i=1}^{i=N} \mathbf{F}_i = N(\mathbf{v}_0 - \mathbf{v}_G) + \nabla \mathbf{v} \cdot \left( \sum_{i=1}^{i=N} L_i \mathbf{p}_i \right) - \left( \sum_{i=1}^{i=N} L_i \dot{\mathbf{p}}_i \right). \quad (35) $$

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

$$ \mathbf{v}_0 = \mathbf{v}_G, \quad (36) $$

implying that the cluster center of gravity is moving with the fluid velocity at that position.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Rigid cluster composed of rods.</figcaption></figure>
The torque created by forces applied on bead i is given by

$$M_i = L_i \mathbf{p}_i \times \mathbf{F}_i. \quad (37)$$

Neglecting inertial effects, the resulting torque for the whole cluster must vanish:

$$\sum_{i=1}^{i=N} \mathbf{M}_i = \mathbf{0}. \quad (38)$$

Taking Eqs. (37) and (34) into account, we have

$$\sum_{i=1}^{i=N} L_i^2 \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i) = \sum_{i=1}^{i=N} L_i^2 \mathbf{p}_i \times \dot{\mathbf{p}}_i. \quad (39)$$

If we define the cluster angular velocity $\omega$ such that

$$\dot{\mathbf{p}}_i = \omega \times \mathbf{p}_i, \quad (40)$$

torque equilibrium reads

$$\sum_{i=1}^{i=N} L_i^2 \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i) = \sum_{i=1}^{i=N} L_i^2 \mathbf{p}_i \times (\omega \times \mathbf{p}_i). \quad (41)$$

Now, by using the vector triple product relationship $\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = \mathbf{b} (\mathbf{a} \cdot \mathbf{c}) - \mathbf{c} (\mathbf{a} \cdot \mathbf{b})$, the right-hand side reads

$$\sum_{i=1}^{i=N} L_i^2 \mathbf{p}_i \times (\omega \times \mathbf{p}_i) = \sum_{i=1}^{i=N} L_i^2 (\omega (\mathbf{p}_i \cdot \mathbf{p}_i) - \mathbf{p}_i (\mathbf{p}_i \cdot \omega)), \quad (42)$$

which, taking into account the normality of vectors $\mathbf{p}_i$ and the fact that $\mathbf{p}_i(\mathbf{p}_i \cdot \omega) = (\mathbf{p}_i \otimes \mathbf{p}_i) \cdot \omega$, becomes

$$\sum_{i=1}^{i=N} L_i^2 \mathbf{p}_i \times (\omega \times \mathbf{p}_i) = \sum_{i=1}^{i=N} L_i^2 (\mathbf{I} - (\mathbf{p}_i \otimes \mathbf{p}_i)) \cdot \omega. \quad (43)$$

with $\mathbf{I}$ the unit matrix.

Now, the left-hand side of Eq. (41) can be rewritten by using the third-order Levi-Civita permutation tensor $\epsilon$ such that $(\mathbf{u} \times \mathbf{v}) = \epsilon : (\mathbf{v} \otimes \mathbf{u})$. We obtain

$$\sum_{i=1}^{i=N} L_i^2 \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i) = \sum_{i=1}^{i=N} L_i^2 \epsilon : (\nabla \mathbf{v} \cdot (\mathbf{p}_i \otimes \mathbf{p}_i)), \quad (44)$$

which finally yields

$$\sum_{i=1}^{i=N} L_i^2 \epsilon : (\nabla \mathbf{v} \cdot (\mathbf{p}_i \otimes \mathbf{p}_i)) = \sum_{i=1}^{i=N} L_i^2 (\mathbf{I} - (\mathbf{p}_i \otimes \mathbf{p}_i)) \cdot \omega. \quad (45)$$

Expressed in a more compact form, we have

$$\epsilon : \nabla \mathbf{v} \cdot \left( \sum_{i=1}^{i=N} L_i^2 (\mathbf{p}_i \otimes \mathbf{p}_i) \right) = \left( \sum_{i=1}^{i=N} L_i^2 \mathbf{I} - \sum_{i=1}^{i=N} L_i^2 (\mathbf{p}_i \otimes \mathbf{p}_i) \right) \cdot \omega. \quad (46)$$

## 3.2. Continuous description

The continuous description considers the pdf $\Upsilon(\mathbf{p}, L)$ giving the fraction of rod beads located at position $\mathbf{p}L$, that can be expressed as

$$\Upsilon(\mathbf{p}, L) = \psi(\mathbf{p}; L) \Gamma(L), \quad (47)$$

where $\psi(\mathbf{p}; L)$ is the angular distribution of beads located at distance $L$ from the cluster center of gravity $\mathcal{G}$, and $\Gamma(L)$ is the fraction of beads located at that distance $L$. The normality condition reads

$$\int_S \psi(\mathbf{p}; L) \, d\mathbf{p} = 1, \quad \forall L. \quad (48)$$

By defining the conformation tensor $\mathbf{c}_L$ related to the population of beads located at distance $L$ with respect to the cluster center of gravity $\mathcal{G}$ as

$$\mathbf{c}_L = \int_S \mathbf{p} \otimes \mathbf{p} \, \psi(\mathbf{p}; L) \, d\mathbf{p}, \quad (49)$$

all sums in the previous expressions must be substituted by the corresponding integrals in the length and orientation domains, $\mathcal{L}$ and $\mathcal{S}$ respectively, weighted by the distribution function $\psi(\mathbf{p}; L) \Gamma(L)$. In particular Eq. (46) becomes:

$$\epsilon : (\nabla \mathbf{v} \cdot \left( \int_{\mathcal{L}} L^2 \left( \int_{S} (\mathbf{p} \otimes \mathbf{p}) \psi(\mathbf{p}; L) \, d\mathbf{p} \right) \Gamma(L) \, dL \right)) = \left( \left( \int_{\mathcal{L}} L^2 \, \Gamma(L) \, dL \right) \mathbf{I} - \left( \int_{\mathcal{L}} L^2 \left( \int_{S} (\mathbf{p} \otimes \mathbf{p}) \psi(\mathbf{p}; L) \, d\mathbf{p} \right) \Gamma(L) \, dL \right) \right) \cdot \omega, \quad (50)$$

which gives, using definition (49),

$$\epsilon : (\nabla \mathbf{v} \cdot \left( \int_{\mathcal{L}} L^2 \mathbf{c}_L \, \Gamma(L) \, dL \right)) = \left( \left( \int_{\mathcal{L}} L^2 \, \Gamma(L) \, dL \right) \mathbf{I} - \left( \int_{\mathcal{L}} L^2 \, \mathbf{c}_L \, \Gamma(L) \, dL \right) \right) \cdot \omega. \quad (51)$$

By introducing the mean square length $\beta$ and tensor $\tilde{\mathbf{c}}$ according to

$$\beta = \int_{\mathcal{L}} L^2 \, \Gamma(L) \, dL, \quad (52)$$

and

$$\tilde{\mathbf{c}} = \int_{\mathcal{L}} L^2 \mathbf{c}_L \, \Gamma(L) \, dL, \quad (53)$$

and due to the fact that $\mathrm{tr}(\mathbf{c}_L) = 1$, we can infer that $\beta$ is in fact the trace of $\tilde{\mathbf{c}}$, that is $\beta = \mathrm{tr}(\tilde{\mathbf{c}})$.

Thus, the particle kinematics finally results:

$$\epsilon : (\nabla \mathbf{v} \cdot \tilde{\mathbf{c}}) = (\mathrm{tr}(\tilde{\mathbf{c}})\mathbf{I} - \tilde{\mathbf{c}}) \cdot \omega, \quad (54)$$

or dividing by $\mathrm{tr}(\tilde{\mathbf{c}})$,

$$\omega = (\mathbf{I} - \mathbf{c})^{-1} \cdot (\epsilon : (\nabla \mathbf{v} \cdot \mathbf{c})), \quad (55)$$

where the conformation tensor $\mathbf{c}$ is given by

$$\mathbf{c} = \frac{\tilde{\mathbf{c}}}{\mathrm{tr}(\tilde{\mathbf{c}})}. \quad (56)$$

**Remark 4.** In the 2D case, $\mathbf{c} \cdot \omega = \mathbf{0}$ and then $\omega = (\epsilon : (\nabla \mathbf{v} \cdot \mathbf{c}))$. 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 $\mathbf{c}$ rotate at the same angular velocity.

As the conformation tensor $\mathbf{c}$ is symmetric and positive definite, it has real eigenvalues and eigenvectors. In 3D, the three mutually perpendicular eigenvectors will be denoted by $\mathbf{u}_1, \mathbf{u}_2$ and $\mathbf{u}_3$, with the associated eigenvalues $\tau_1, \tau_2$ and $\tau_3$ respectively.

Imagine a rigid cluster composed of three rods oriented in directions $\mathbf{u}_1, \mathbf{u}_2$ and $\mathbf{u}_3$ with respective lengths $\sqrt{\tau_1}, \sqrt{\tau_2}$ and $\sqrt{\tau_3}$. The conformation tensor of such a three-rod cluster coincides with $\mathbf{c}$ and then both tensors have the same rotary velocity.

If two eigenvalues are identical, e.g. $\tau_2 = \tau_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 = \sqrt{\frac{\tau_1}{\tau_2}}$ given by Jeffery's Eq. (2).
# 4. Numerical results

## 4.1. 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 $\infty$: $r \to \infty$ (rods), $r > 1$ (prolates), $r = 1$ (sphere), $0 < r < 1$ (oblates) and $r \to 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.

<caption>**Table 1**</caption>

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

<table><thead><tr><th>N</th><th>M</th><th>r</th><th>r<sub>e</sub></th><th>r<sub>e</sub><sup>zh</sup></th></tr></thead><tbody><tr><td>4</td><td>3</td><td>1</td><td>1.15</td><td>1.04</td></tr><tr><td>4</td><td>4</td><td>1</td><td>1.05</td><td>1.04</td></tr><tr><td>4</td><td>5</td><td>1</td><td>1.00</td><td>1.04</td></tr><tr><td>4</td><td>6</td><td>1</td><td>0.97</td><td>1.04</td></tr><tr><td>6</td><td>3</td><td>1</td><td>1.15</td><td>1.04</td></tr><tr><td>8</td><td>3</td><td>1</td><td>1.15</td><td>1.04</td></tr><tr><td>4</td><td>3</td><td>5</td><td>5.77</td><td>4.53</td></tr><tr><td>4</td><td>6</td><td>5</td><td>4.83</td><td>4.53</td></tr><tr><td>4</td><td>3</td><td>10</td><td>11.54</td><td>8.66</td></tr></tbody></table>

## 4.2. 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.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Angular velocity distribution for $r = 1$ and $M = 3, 4, 6$.</figcaption></figure>
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:

$$r_e^{zh} = 0.000035 \cdot r^3 - 0.00467 \cdot r^2 + 0.764 \cdot r + 0.404. \quad (57)$$

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 \cdot M$ beads $B_{ij}$, with coordinates $x_i = \frac{D}{2} \cdot \cos(\varphi_i)$, $y_i = \frac{D}{2} \cdot \sin(\varphi_i)$ and $z_j = -\frac{H}{2} + h_z \cdot (j-1)$ ($i = 1, \dots, N$ and $j = 1, \dots, M$), with $\varphi_i = h_\varphi \cdot (i-1)$, $h_\varphi = \frac{2\pi}{N}$ and $h_z = \frac{H}{M-1}$.

We study different configurations by selecting different values of $N$ and $M$, computing the associated conformation tensors $\mathbf{c}$ from

$$\mathbf{c} = \frac{\sum_{i=1}^{i=N} L_i^2 (\mathbf{p}_i \otimes \mathbf{p}_i)}{\sum_{i=1}^{i=N} L_i^2}, \quad (58)$$

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).

Fig. 3 illustrates the case of a simple shear flow $\mathbf{v}^T = (0, \dot{\gamma}, 0)$ with $\dot{\gamma} = 1s^{-1}$. We show the angular velocity $\omega_x$ for different cylin-

<figure><img src="image_5.png" /><figcaption>Fig. 5. Angular velocity distribution for $r = 10$ and $M = 3$.</figcaption></figure>

der orientations $\mathbf{p}(\theta) = (0, \cos(\theta), \sin(\theta))$ for $r = 1, N = 4$ and $M = 3, 4, 6$.

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 $\lambda$ 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 \to \infty$).

## 5. 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 $\mathbf{c}$ with unit trace that can be considered as the appropriate conformation tensor for describing cluster kinematics. When the conformation tensor $\mathbf{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

<figure><img src="image_6.png" /><figcaption>Fig. 4. Angular velocity distribution for $r = 5$ and $M = 3, 6$.</figcaption></figure>
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 \approx 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 $\omega$, we can obtain the expression of $\dot{\mathbf{p}}$ in order to define the conformation evolution from $\dot{\mathbf{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 $\mathbf{c}_i$, $i = 1, \dots, N$. Within the continuous framework, the suspension is characterized by the pdf $\Psi$, which now depends on the physical coordinates (space **x** and time **t**) and the conformation coordinate **c**. Thus, the pdf reads $\Psi(\mathbf{x}, t, \mathbf{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 $\dot{\mathbf{c}}_i$. Within the continuous framework, the pdf $\Psi(\mathbf{x}, t, \mathbf{c})$ evolves according to its associated Fokker-Planck equation, with the knowledge of $\dot{\mathbf{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.

7. The macroscopic description uses a coarser description based on the moments of the pdf. Here, the simplest choice consists of the first moment $\mathbf{C}(\mathbf{x}, t)$ defined as

$$ \mathbf{C}(\mathbf{x}, t) = \int_C \mathbf{c} \, \Psi(\mathbf{x}, t, \mathbf{c}) \, d\mathbf{c}. \qquad (59) $$

8. 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.

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 $\mathbf{m} = \mathbf{I} - \mathbf{c}$, as matrix **m** is $3 \times 3$ it satisfies its characteristic polynomial (according to the Cayley-Hamilton theorem):

$$ \mathbf{m}^3 - tr(\mathbf{m})\mathbf{m}^2 + \frac{(tr(\mathbf{m}))^2 - tr(\mathbf{m}^2)}{2} \mathbf{m} - det(\mathbf{m})\mathbf{I} = \mathbf{0}, \qquad (A.1) $$

where $det(\mathbf{m})$ refers to the matrix determinant.

Now we can write Eq. (A.1) in the form

$$ \mathbf{m} \ \mathbf{m}^2 - tr(\mathbf{m})\mathbf{m} + \frac{(tr(\mathbf{m}))^2 - tr(\mathbf{m}^2)}{2} \mathbf{I} = det(\mathbf{m})\mathbf{I}, \qquad (A.2) $$

from which we can infer that $\mathbf{m}^{-1}$ is given by

$$ \mathbf{m}^{-1} = \frac{1}{det(\mathbf{m})} \ \mathbf{m}^2 - tr(\mathbf{m})\mathbf{m} + \frac{(tr(\mathbf{m}))^2 - tr(\mathbf{m}^2)}{2} \mathbf{I}. \qquad (A.3) $$

By developing the expression of $\mathbf{m}^{-1}$ and taking into account that $\mathbf{m} = \mathbf{I} - \mathbf{c}$, we obtain

$$ (\mathbf{I} - \mathbf{c})^{-1} = \alpha_1 \mathbf{I} + \alpha_2 \mathbf{c}^2, \qquad (A.4) $$

where the coefficients $\alpha_1$ and $\alpha_2$, that depend on the invariants of **c**, are given by

$$ \alpha_1 = \frac{\frac{1}{2}(\mathbf{I} - \mathbf{c} : \mathbf{c})}{\frac{1}{2}(\mathbf{I} - \mathbf{c} : \mathbf{c}) - det(\mathbf{c})}, \qquad (A.5) $$

and

$$ \alpha_2 = \frac{1}{\frac{1}{2}(\mathbf{I} - \mathbf{c} : \mathbf{c}) - det(\mathbf{c})}. \qquad (A.6) $$

Thus, finally Eq. (55) reads

$$ \omega = (\alpha_1 \mathbf{I} + \alpha_2 \mathbf{c}^2) \cdot (\epsilon : (\nabla \mathbf{v} \cdot \mathbf{c})). \qquad (A.7) $$

# Appendix B. Jeffery revisited

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]:

$$ \dot{\mathbf{p}} = \mathbf{\Omega} \cdot \mathbf{p} + \lambda \cdot \mathbf{D} \cdot \mathbf{p} - \lambda (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p}. \qquad (B.1) $$

where $\lambda$ depends on the ellipsoid aspect ratio $r$ (the length to diameter ratio):

$$ \lambda = \frac{r^2 - 1}{r^2 + 1}, \qquad (B.2) $$

and where for the sake of clarity we write the term **D**: $(\mathbf{p} \otimes \mathbf{p})$ in the matrix algebra format, by using the equivalence **D**: $(\mathbf{p} \otimes \mathbf{p}) \equiv \mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{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 \ge L_2$, as depicted in Fig. 6, such that the parameter $\lambda$ in Jeffery's Eq. (B.1) results from $r = \frac{L_1}{L_2}$,

$$ \lambda = \frac{L_1^2 - L_2^2}{L_1^2 + L_2^2}. \qquad (B.3) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Two-rod rigid cluster representing an ellipsoidal 2D particle.</figcaption></figure>

In the present configuration the forces applied on bead $L_1\mathbf{p}_1$ and $L_2\mathbf{p}_2$ are

$$ \mathbf{F}_1 = \xi L_1 (\nabla \mathbf{v} \cdot \mathbf{p}_1 - \dot{\mathbf{p}}_1), \quad (B.4) $$

and

$$ \mathbf{F}_2 = \xi L_2 (\nabla \mathbf{v} \cdot \mathbf{p}_2 - \dot{\mathbf{p}}_2), \quad (B.5) $$

respectively, with

$$ \begin{cases} \dot{\mathbf{p}}_1 = \omega \times \mathbf{p}_1 \\ \dot{\mathbf{p}}_2 = \omega \times \mathbf{p}_2. \end{cases} \quad (B.6) $$

Torque balance implies:

$$ L_1^2 \mathbf{p}_1 \times (\nabla \mathbf{v} \cdot \mathbf{p}_1 - \dot{\mathbf{p}}_1) + L_2^2 \mathbf{p}_2 \times (\nabla \mathbf{v} \cdot \mathbf{p}_2 - \dot{\mathbf{p}}_2) = \mathbf{0}. \quad (B.7) $$

Introducing (B.6) into (B.7) and taking into account that $\mathbf{p}_1 \times \omega \times \mathbf{p}_1 = \omega$ and $\mathbf{p}_2 \times \omega \times \mathbf{p}_2 = \omega$, we obtain

$$ \omega = \frac{L_1^2}{L_1^2 + L_2^2} (\mathbf{p}_1 \times (\nabla \mathbf{v} \cdot \mathbf{p}_1)) + \frac{L_2^2}{L_1^2 + L_2^2} (\mathbf{p}_2 \times (\nabla \mathbf{v} \cdot \mathbf{p}_2)). \quad (B.8) $$

Now, the ellipsoid rotary velocity is given by

$$ \begin{aligned} \dot{\mathbf{p}}_1 &= \omega \times \mathbf{p}_1 \\ &= \frac{L_1^2}{L_1^2 + L_2^2} ((\mathbf{p}_1 \times (\nabla \mathbf{v} \cdot \mathbf{p}_1)) \times \mathbf{p}_1) \\ &\quad + \frac{L_2^2}{L_1^2 + L_2^2} ((\mathbf{p}_2 \times (\nabla \mathbf{v} \cdot \mathbf{p}_2)) \times \mathbf{p}_1). \end{aligned} \quad (B.9) $$

It is easy to verify that $\dot{\mathbf{p}}_1 \cdot \mathbf{p}_1 = 0$ as expected, i.e. the variation of $\mathbf{p}_1$ is directed along its perpendicular direction, which in this case corresponds to $\mathbf{p}_2$.

Now, applying again the triple vector product formula $(\mathbf{a} \times \mathbf{b}) \times \mathbf{c} = -\mathbf{a} \cdot (\mathbf{b} \cdot \mathbf{c}) + \mathbf{b} \cdot (\mathbf{a} \cdot \mathbf{c})$, Eq. (B.9) reads

$$ \begin{aligned} \dot{\mathbf{p}}_1 = & \frac{L_1^2}{L_1^2 + L_2^2} (\nabla \mathbf{p}_1 - (\mathbf{p}_1^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_1) \mathbf{p}_1) \\ & - \frac{L_1^2}{L_1^2 + L_2^2} ((\mathbf{p}_1^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_2) \mathbf{p}_2). \end{aligned} \quad (B.10) $$

We develop the last term of Eq. (B.10) to obtain an equation that only contains $\mathbf{p}_1$ in order to compare with Jeffery's expression (B.1).

First, we apply the decomposition $\nabla \mathbf{v} = \mathbf{D} + \mathbf{\Omega}$ from which the last term in Eq. (B.10) reads

$$ (\mathbf{p}_1^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_2) \mathbf{p}_2 = (\mathbf{p}_1^T \cdot \mathbf{D} \cdot \mathbf{p}_2) \mathbf{p}_2 + (\mathbf{p}_1^T \cdot \mathbf{\Omega} \cdot \mathbf{p}_2) \mathbf{p}_2. \quad (B.11) $$

We proceed from the first term in the right-hand side of Eq. (B.11), using the fact that $\mathbf{D}$ is symmetric:

$$ \begin{aligned} (\mathbf{p}_1^T \cdot \mathbf{D} \cdot \mathbf{p}_2) \mathbf{p}_2 &= \mathbf{p}_2 (\mathbf{p}_1^T \cdot \mathbf{D} \cdot \mathbf{p}_2) = \mathbf{p}_2 (\mathbf{p}_2^T \cdot \mathbf{D} \cdot \mathbf{p}_1) \\ &= (\mathbf{p}_2 \otimes \mathbf{p}_2) \cdot \mathbf{D} \cdot \mathbf{p}_1. \end{aligned} \quad (B.12) $$

Since $\mathbf{p}_1$ and $\mathbf{p}_2$ are mutually perpendicular, we have

$$ (\mathbf{p}_2 \otimes \mathbf{p}_2) + (\mathbf{p}_1 \otimes \mathbf{p}_1) = \mathbf{I}, \quad (B.13) $$

and Eq. (B.12) can then be written as

$$ (\mathbf{p}_1^T \cdot \mathbf{D} \cdot \mathbf{p}_2) \mathbf{p}_2 = (\mathbf{I} - \mathbf{p}_1 \otimes \mathbf{p}_1) \cdot \mathbf{D} \cdot \mathbf{p}_1. \quad (B.14) $$

Finally, we consider the second term in the right-hand side of Eq. (B.11), using the fact that $\mathbf{\Omega}$ is skew-symmetric:

$$ \begin{aligned} (\mathbf{p}_1^T \cdot \mathbf{\Omega} \cdot \mathbf{p}_2) \mathbf{p}_2 &= \mathbf{p}_2 (\mathbf{p}_1^T \cdot \mathbf{\Omega} \cdot \mathbf{p}_2) = -\mathbf{p}_2 (\mathbf{p}_2^T \cdot \mathbf{\Omega} \cdot \mathbf{p}_1) \\ &= -(\mathbf{p}_2 \otimes \mathbf{p}_2) \cdot \mathbf{\Omega} \cdot \mathbf{p}_1, \end{aligned} \quad (B.15) $$

which finally becomes

$$ (\mathbf{p}_1^T \cdot \mathbf{\Omega} \cdot \mathbf{p}_2) \mathbf{p}_2 = -(\mathbf{I} - \mathbf{p}_1 \otimes \mathbf{p}_1) \cdot \mathbf{\Omega} \cdot \mathbf{p}_1. \quad (B.16) $$

Now, coming back to Eq. (B.10), we obtain

$$ \dot{\mathbf{p}}_1 = \mathbf{\Omega} \cdot \mathbf{p}_1 + \frac{L_1^2 - L_2^2}{L_1^2 + L_2^2} \mathbf{D} \cdot \mathbf{p}_1 - \frac{L_1^2 - L_2^2}{L_1^2 + L_2^2} (\mathbf{p}_1^T \cdot \mathbf{D} \cdot \mathbf{p}_1) \mathbf{p}_1, \quad (B.17) $$

which corresponds exactly to Jeffery's expression for a 2D ellipsoid of aspect ratio $L_1/L_2$.

For circular particles characterized by a unit aspect ratio $L_1 = L_2$, we obtain as expected

$$ \dot{\mathbf{p}} = \mathbf{\Omega} \cdot \mathbf{p}. \quad (B.18) $$

## References

[1] E. Abisset-Chavanne, R. Mezher, S. Le Corre, A. Ammar, F. Chinesta, Kinetic theory microstructure modeling in concentrated suspensions, Entropy 15 (2013) 2805-2832.

[2] S. Advani, Ch. Tucker, The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. 31 (1987) 751-784.

[3] S. Advani, Ch. Tucker, Closure approximations for three-dimensional structure tensors, J. Rheol. 34 (1990) 367-386.

[4] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, J. Non-Newtonian Fluid Mech. 139 (2006) 153-176.

[5] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: Transient simulation using space-time separated representations, J. Non-Newtonian Fluid Mech. 144 (2007) 98-121.

[6] G.K. Batchelor, The stress system in a suspension of force-free particles, J. Fluid Mech. 41 (1970) 545-570.

[7] R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamic of Polymeric Liquid, Kinetic Theory, vol. 2, John Wiley and Sons, 1987.

[8] H. Brenner, Rheology of a dilute suspension of axisymmetric Brownian particles, Int. J. Multiphase Flow 1 (1974) 195-341.

[9] F.P. Bretherton, The motion of rigid particles in a shear flow at low Reynolds number, J. Fluid Mech. 14 (1962) 284-304.

[10] J.M. Burgers, On the motion of small particles of elongated form suspended in a viscous liquid, in: Second Report on Viscosity and Plasticity, Kon. Ned. Akad. Wet., Verhand. (Erste Sectie), North-Holland Publ., Amsterdam, 1938, pp. 113-184 (Chapter 3).

[11] C.V. Chaubal, A. Srinivasan, O. Egecioglu, L.G. Leal, Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems, J. Non-Newtonian Fluid Mech. 70 (1997) 125-154.

[12] C. Chauviere, A. Lozinski, Simulation of dilute polymer solutions using a Fokker-Planck equation, Comput. Fluids 33 (2004) 687-696.

[13] K. Chiba, A. Ammar, F. Chinesta, On the fiber orientation in steady recirculating flows involving short fibers suspensions, Rheol. Acta 44 (2005) 406-417.

[14] F. Chinesta, G. Chaidron, A. Poitou, On the solution of the Fokker-Planck equation in steady recirculating flows involving short fibre suspensions, J. Non-Newtonian Fluid Mech. 113 (2003) 97-125.

[15] F. Chinesta, A. Ammar, A. Leygue, R. Keunings, An overview of the proper generalized decomposition with applications in computational rheology, J. Non Newtonian Fluid Mech. 166 (2011) 578-592.
[16] F. Chinesta, From single-scale to two-scales kinetic theory descriptions of rods suspensions, Arch. Comput. Methods Eng. 20 (1) (2013) 1–29.

[17] F. Chinesta, R. Keunings, A. Leygue, The Proper Generalized Decomposition for Advanced Numerical Simulations, A Primer. Springerbriefs, Springer, 2014.

[18] R.G. Cox, The motion of long slender bodies in a viscous fluid. Part 2. Shear flow, J. Fluid Mech. 45 (1971) 625–657.

[19] C. Cruz, L. Illoul, F. Chinesta, G. Regnier, Effects of a bent structure on the linear viscoelastic response of carbon nanotube diluted suspensions, Rheol. Acta 49 (2010) 1141–1155.

[20] C. Cruz, F. Chinesta, G. Regnier, Review on the Brownian dynamics simulation of bead-rod-spring models encountered in computational rheology, Arch. Comput. Methods Eng. 19 (2) (2012) 227–259.

[21] M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Clarendon Press, Oxford, 1987.

[22] F. Dupret, V. Verleye, Modelling the flow of fibre suspensions in narrow gaps, in: D.A. Siginer, D. De Kee, R.P. Chabra (Eds.), Advances in the Flow and Rheology of Non-Newtonian Fluids, Rheology Series, Elsevier, 1999, pp. 1347–1398.

[23] F. Folgar, Ch. Tucker, Orientation behavior of fibers in concentrated suspensions, J. Reinf. Plast. Comp. 3 (1984) 98–119.

[24] G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech. 13 (1962) 33–62.

[25] J. Hinch, G. Leal, The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles, J. Fluid Mech. 52 (1972) 683–712.

[26] J. Hinch, G. Leal, Constitutive equations in suspension mechanics. Part I, J. Fluid Mech. 71 (1975) 481–495.

[27] J. Hinch, G. Leal, Constitutive equations in suspension mechanics. Part II, J. Fluid Mech. 76 (1976) 187–208.

[28] J. Hinch, G. Leal, Rotation of small non-axisymmetric particles in a simple shear flow, J. Fluid Mech. 92 (1979) 591–608.

[29] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proc. R. Soc. Lond. A 102 (1922) 161–179.

[30] R. Keunings, On the Peterlin approximation for finitely extensible dumbbells, J. Non-Newtonian Fluid Mech. 68 (1997) 85–100.

[31] R. Keunings, Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory, in: D.M. Binding, K. Walters (Eds.), Rheology Reviews, British Society of Rheology, 2004, pp. 67–98.

[32] M. Kroger, A. Ammar, F. Chinesta, Consistent closure schemes for statistical models of anisotropic fluids, J. Non-Newtonian Fluid Mech. 149 (2008) 40–55.

[33] A. Lozinski, C. Chauviere, A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model, J. Comput. Phys. 189 (2003) 607–625.

[34] A. Ma, F. Chinesta, M. Mackley, A. Ammar, The rheological modelling of Carbon Nanotube (CNT) suspensions in steady shear flows, Int. J. Mater. Form. 2 (2008) 83–88.

[35] A. Ma, F. Chinesta, A. Ammar, M. Mackley, Rheological modelling of carbon nanotube aggregate suspensions, J. Rheol. 52 (6) (2008) 1311–1330.

[36] A. Ma, F. Chinesta, M. Mackley, The rheology and modelling of chemically treated carbon nanotube suspensions, J. Rheol. 53 (3) (2009) 547–573.

[37] B. Mokdad, E. Pruliere, A. Ammar, F. Chinesta, On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach, Appl. Rheol. 17 (2) (2007) 1–14. 26494.

[38] H.C. Öttinger, M. Laso, Smart polymers in finite element calculation, Int. Congr. on Rheology, Brussel, Belguim, 1992.

[39] C. Petrie, The rheology of fibre suspensions, J. Non-Newtonian Fluid Mech. 87 (1999) 369–402.

[40] E. Pruliere, A. Ammar, N. El Kissi, F. Chinesta, Recirculating flows involving short fiber suspensions: numerical difficulties and efficient advanced micro-macro solvers, Arch. Comput. Methods Eng., State Art Rev. 16 (2009) 1–30.

[41] B.J. Trevelyan, S.G. Mason, Particle motions in sheared suspensions. I. Rotations, J. Colloid Sci. 6 (1951) 354–367.

[42] P. Wapperom, R. Keunings, Numerical simulation of branched polymer melts in transient complex flows using pom-pom models, J. Non-Newtonian Fluid Mech. 97 (2001) 267–281.

[43] D. Zhang, D.E. Smith, D.A. Jack, S. Montgomery-Smith, Numerical evaluation of single fiber motion for short-fiber-reinforced composite materials processing, J. Manuf. Sci. Eng. 133 (5) (2011) 051002.