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

HAL
open science

# From Single-Scale to Two-Scales Kinetic Theory Descriptions of Rods Suspensions

Francisco Chinesta

► To cite this version:

Francisco Chinesta. From Single-Scale to Two-Scales Kinetic Theory Descriptions of Rods Suspensions. Archives of Computational Methods in Engineering, 2013, 20 (1), pp.1-29. 10.1007/s11831-013-9079-3. http://dx.doi.org/10.1007/s11831-013-9079-3

HAL Id: hal-01006908

https://hal.science/hal-01006908

Submitted on 29 Apr 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.

Public Domain
# From single-scale to two-scales kinetic theory descriptions of rods suspensions

Francisco Chinesta

**Abstract** This paper proposes a first attempt to define a two-scales kinetic theory description of suspensions involving short fibers, nano-fibers or nanotubes. We start revisiting the description of dilute enough suspensions for which microscopic, mesoscopic and macroscopic descriptions are available and all them have been successfully applied for describing the rheology of such suspensions. When the suspensions become too concentrated fiber-fiber interactions cannot be neglected and then classical dilute theories fail for describing the rich microstructure evolution. In the semi-concentrated regime some interaction mechanisms that mimetic the randomizing effect of fiber-fiber interactions were successfully introduced. Finally, when the concentration becomes high enough, richer microstructures can be observed. They involve a diversity of fiber clusters or aggregates with complex kinematics, and different sizes and shapes. These clusters can interact to create larger clusters and also break because the flow induced hydrodynamic forces. In this paper we propose a double-scale kinetic theory model that at the first scale consider the kinematics of the clusters, whose structure itself is described at the finest scale, the one related to the rods constituting the clusters.

## 1 Introduction

Short fibers, nanofibers or CNTs suspensions present different morphologies depending on their concentrations. When the concentration is dilute enough one can describe the microstructure by tracking a population of rods that move with the suspending fluid and orient depending on the velocity gradient according to the Jeffery's equation [38] that relates the orientation evolution with the flow velocity field. In that case the motion and orientation of each fiber is assumed decoupled from the others. The representative

This work has been partially supported by the IUF – *Institut Universitaire de France* –.

Francisco Chinesta
EADS Foundation Chair "Advanced Computational Manufacturing Processes"
GEM, UMR CNRS - Centrale Nantes
Institut Universitaire de France
1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France
E-mail: Francisco.Chinesta@ec-nantes.fr
population can involve too many fibers, and in that case the computational efforts to track the population is unaffordable. Thus, the simple and well-defined physics must be sacrificed in order to derive coarser descriptions.

Kinetic theory approaches [27] [15] describe such systems at the mesoscopic scale. Their main advantage is their capability to address macroscopic systems, while keeping the fine physics through a number of conformational coordinates introduced for describing the microstructure and its time evolution. At this mesoscopic scale, the microstructure is defined from a distribution function that depends on the physical space, the time and a number of conformational coordinates - the rods orientation in the case of slender bodies suspensions -. The moments of this distribution constitutes a coarser description in general used in macroscopic modeling [2].

At the macroscopic scale the equations governing the time evolution of these moments usually involve closure approximations whose impact in the results is unpredictable. Other times, macroscopic equations are carefully postulated in order to guarantee the model objectivity and the thermodynamical admissibility.

In the case of dilute suspensions of short fibers the three scales have been extensively considered to model such systems, without major difficulties. However, as soon as the concentration increases the difficulties appear. In the semi-dilute or semi-concentrated regimes fiber-fiber interactions occur, but in general they can be accurately modeled by introducing a sort of randomizing diffusion term. Again, models based on this mechanism describe quite well the experimental observations. The reader can refer to [4] [51] [34] [13] [35] [36] [37] [54] and the references therein. The worst scenario concerns concentrated suspensions involving entangled clusters of particles and in which aggregation/disaggregation takes place.

The first natural question is how describing such systems? At the macroscopic scale one could try to fit some power-law constitutive equation, however, this description does not allow to describe the microstructure. At the microscopic scale, direct numerical simulations describing complex fiber-fiber interactions can be carried out in small enough representative volumes. A natural candidate to be a reasonable compromise between (fine) micro and (fast) macro descriptions consists of considering again a kinetic theory description.

The main issue of such an approach lies in the fact that it must include two scales, the one involving the aggregates and the one related to the rods constituting the aggregates. What are the appropriate conformational coordinates? How to determine the time evolution of these conformational coordinates? How to represent simultaneously both physical scales, one related to the aggregates and the other to the fibers? How to derive the aggregation/disaggregation mechanisms? How to derive the expression of the stress tensor? ... Some key questions that we will address in this work.

## 1.1 A concrete example: suspensions of carbon nanotubes

All these questions arise for example when addressing carbon nanotubes suspensions the more and more used in nanocomposites [45] [26].

In the case of chemically treated single-walled carbon nanotubes (CNTs) suspended within an epoxy resin the aggregation is prevented and we can consider that we are dealing with a large population of free rods in the dilute regime or rods experiencing interactions in the semi-concentrated or concentrated regimes. Thus, when a CNT suspension was subject to a steady shear flow, it exhibited a shear-thinning characteristic,
which was subsequently modelled by a Fokker-Planck (FP) based orientation model
[47]. The model assumes that shear flow aligns CNT in the flow direction, but there
are events such as Brownian motion and tube-tube interaction trying to randomise the
orientation. In the FP orientation model, randomising events were modelled with an
appropriate rotary diffusion coefficient *D*<sub>r</sub> and the shear-thinning behaviour was ex-
plained in terms of progressive alignment of CNTs towards the shear direction. In terms
of linear viscoelasticity (LVE), small-amplitude oscillatory measurements revealed mild
elasticity for semi-dilute treated CNT suspensions. The exact origin for this elasticity
is not clear and both tube-tube interaction and bending/stretching of CNTs have been
proposed by other authors as possible origins (see [24] and the references therein).

In [47] experimental LVE data of the treated CNT suspensions were fitted using
the FP orientation model with an "effective diffusion coefficient" term and an empirical
relation was subsequently identified for the effective diffusion term. Intuitively, chemical
treatment has created a weakly interconnected network of CNT and it is believed that
the mild elasticity originated from this weak network as well as other randomizing
events (Brownian motion and tube-tube hydrodynamic interaction). Finally, step strain
experiments confirmed the presence of a weak network at small strains, which at large
strains was found to be destroyed. Incorporation of a strain softening factor allowed
for the formulation of a self-consistent FP based orientation model describing both
the steady shear and LVE responses of treated CNT suspensions. A finer Brownian
dynamics modeling was addressed in [25], where the elasticity effects were explained as
a direct consequence of the nanotube bending, more in particular, the bending of those
having a non-straight natural configuration because the presence of side-wall defects
at the origin of the CNT bent conformation.

The rheological modelling of untreated Carbon Nanotubes (CNTs) suspended within
an epoxy resin in a kinetic theory framework was addressed in [?]. The untreated CNT
suspensions exhibited significant steady shear-thinning and contained optically resolv-
able aggregate structures depending on the applied shear rate as observed in Fig. 1.
A simple orientation model, based on a Fokker-Planck advection-diffusion description,
failed to capture the experimentally observed rheological responses for untreated CNT
suspensions. A new model named the "Aggregation/Orienta tion" (AO) model has been
developed to describe the experimental findings. The model integrated elements of both
a standard orientation model and aggregation modelling concepts within the Fokker-
Planck formalism. A hierarchy of states between CNTs that are free from entanglement
and a complete CNT network was incorporated into the AO model, thereby enabling
different microstructure populations to exist for different shear conditions. Using a
small number of adjustable parameters, it was found that the experimental data could
be fitted with reasonable precision. However, important assumptions, somewhat unjus-
tified, concerning the microstructure evolution of CNTs clusters immersed in the flow
were introduced.

In this work we are first revisiting the description of dilute fibers suspensions, before
moving to more concentrated regimes exhibiting richer evolving microstructures.

## 2 Dilute suspensions of non-Brownian rods

We consider a suspending medium consisting of a newtonian fluid of viscosity η in which
there are suspended N rigid slender rods. We assume as first approximation that the
fibers presence and orientation do not affect the flow kinematics that is defined by the
<figure><img src="image_2.png" /><figcaption>Fig. 1 Optical microstructure of 0.1% CNT suspended within an epoxy resin after shearing at (a) 0 s⁻¹, (b) 0.5 s⁻¹, (c) 10 s⁻¹ and (d) 100 s⁻¹. The micrographs were captured after a shearing time of 100 seconds and using the Cambridge Shear System (CSS 450) with an optical depth of 130 µm.</figcaption></figure>

velocity field $\mathbf{v}(\mathbf{x}, t)$, with $\mathbf{x} \in \Omega \in \mathbb{R}^d$. For the sake of simplicity, and without loss of generality, in what follows we are considering a 2D domain, i.e. $\Omega \in \mathbb{R}^2$.

## 2.1 Microscopic description

The microstructure could be described at the microscopic scale by given the unit vector defining the orientation of each rod, i.e. $\mathbf{p}_i$, $i = 1, \dots, N$. For ellipsoids the orientation evolution is given by the Jeffery's equation [38] that particularized for rods (ellipsoids with infinite aspect ratio) leads to:

$$ \dot{\mathbf{p}}_i = \nabla \mathbf{v} \cdot \mathbf{p}_i - (\mathbf{p}_i^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_i) \cdot \mathbf{p}_i, \quad i = 1, \dots, N \quad (1) $$

It can be noticed that this expression consists of two contributions, a sort of affine deformation given by the term $\nabla \mathbf{v} \cdot \mathbf{p}_i$ and the second one that removes the deformation in the rod direction because the rod is considered inextensible. It easy to verify that $\mathbf{p}_i^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_i \equiv \mathbf{p}_i^T \cdot \mathbf{D} \cdot \mathbf{p}_i$, being $\mathbf{D}$ the rate of strain tensor (symmetric component of the gradient of velocity tensor $\nabla \mathbf{v}$).

This equation could be derived by considering the system illustrated in Fig. 2 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 $\mathbf{F}_i$ depend on the difference of velocities between the fluid and the bead, the first one given by $\nabla \mathbf{v} \cdot \mathbf{p}_i \cdot L$ and the second one by $\dot{\mathbf{p}}_i \cdot L$. Thus, forces $\mathbf{F}_i$ write:
<figure><img src="image_3.png" /><figcaption>Fig. 2 Hydrodynamic forces applying on a rod immersed in a Newtonian fluid</figcaption></figure>

$$ \mathbf{F}_i = \xi \cdot (\nabla \mathbf{v} \cdot \mathbf{p}_i \cdot L - \dot{\mathbf{p}}_i \cdot L) \quad (2) $$

being $\xi$ the friction coefficient.

Obviously if $\mathbf{F}_i$ applies on the bead $\mathbf{p}_i \cdot L$, then in the opposite bead $-\mathbf{p}_i \cdot L$ applies $-\mathbf{F}_i$ ensuring a null balance of forces compatible with the assumed rood inertialess. However, because the rod mass is assumed null, and then its associated inertia moment, the resulting torque must also vanish. The only possibility for this is that force $\mathbf{F}_i$ acts along $\mathbf{p}_i$, that is $\mathbf{F}_i = \lambda \cdot \mathbf{p}_i$, with $\lambda \in \mathcal{R}$. Thus we can write

$$ \lambda \cdot \mathbf{p}_i = \xi \cdot (\nabla \mathbf{v} \cdot \mathbf{p}_i \cdot L - \dot{\mathbf{p}}_i \cdot L) \quad (3) $$

Premultiplying Eq. (3) by $\mathbf{p}_i$ and taking into account that $\mathbf{p}_i^T \cdot \mathbf{p}_i = 1$ and consequently $\mathbf{p}_i^T \cdot \dot{\mathbf{p}}_i = 0$, it results:

$$ \lambda = \xi \cdot L \cdot \mathbf{p}_i^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_i \quad (4) $$

implying

$$ \mathbf{F}_i = \xi \cdot L \cdot (\mathbf{p}_i^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_i) \cdot \mathbf{p}_i \quad (5) $$

that introduced into Eq. (2) leads to the Jeffery's equation for infinite aspect ratio ellipsoids:

$$ \dot{\mathbf{p}}_i = \nabla \mathbf{v} \cdot \mathbf{p}_i - (\mathbf{p}_i^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_i) \cdot \mathbf{p}_i \quad (6) $$

The extra-stress tensor results from the addition of the suspending medium contribution $\boldsymbol{\tau}^f$ and the contribution due to the forces applying on the rods $\boldsymbol{\tau}^r$. The first contribution can be expressed from $\boldsymbol{\tau}^f = 2 \cdot \eta \cdot \mathbf{D}$, being the second one

$$ \boldsymbol{\tau}^r = L \cdot \sum_{i=1}^{i=N} \mathbf{p}_i \otimes \mathbf{F}_i = L \cdot \sum_{i=1}^{i=N} \lambda \cdot \mathbf{p}_i \otimes \mathbf{p}_i \quad (7) $$

where the symbol "$\otimes$" denotes the tensor product that for two vectors $\mathbf{a}$ and $\mathbf{b}$ results $(\mathbf{a} \otimes \mathbf{b})_{ij} = \mathbf{a}_i \cdot \mathbf{b}_j$.

Introducing the expression of $\lambda$ it results

$$ \boldsymbol{\tau}^r = \sum_{i=1}^{i=N} \xi \cdot L^2 \cdot (\mathbf{p}_i^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}_i) \cdot (\mathbf{p}_i \otimes \mathbf{p}_i) \quad (8) $$
Eq. (8) can be rewritten as:

$$ \boldsymbol{\tau}^r = \sum_{i=1}^{i=N} \xi \cdot L^2 \cdot ((\nabla \mathbf{v}) : (\mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i)) \quad (9) $$

where the symbol " : " denotes the tensor product twice contracted.

The Cauchy stress tensor reads:

$$ \boldsymbol{\tau} = \boldsymbol{\tau}^f + \boldsymbol{\tau}^r = 2 \cdot \eta \cdot \mathbf{D} + \sum_{i=1}^{i=N} \xi \cdot L^2 \cdot ((\nabla \mathbf{v}) : (\mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i)) \quad (10) $$

It is important to notice that as forces apply in the rod direction, the resulting stress tensor will be symmetric, i.e. $\boldsymbol{\tau} = \boldsymbol{\tau}^T$.

Thus using all the previous developments, the microstructure evolution could be described as follows:

- Given the flow velocity field $\mathbf{v}(\mathbf{x}, t)$, the initial position of the centre of gravity of each rod $\mathbf{x}_i^G(t = 0)$ and its orientation $\mathbf{p}_i(t = 0)$ or knowing both them at time $t$

- For each fiber $i = 1, \dots, N$:

- Update the position of the centre of gravity by integrating:

$$ \frac{d\mathbf{x}_i^G}{dt} = \mathbf{v}(\mathbf{x}_i^G, t) \quad (11) $$

The simplest integration, by considering a time step $\Delta t$, consists in the first order backward scheme that writes:

$$ \mathbf{x}_i^G(t + \Delta t) \approx \mathbf{x}_i^G(t) + \mathbf{v}(\mathbf{x}_i^G, t) \cdot \Delta t \quad (12) $$

- Update its orientation by integrating

$$ \frac{d\mathbf{p}_i}{dt} = \nabla\mathbf{v}(\mathbf{x}_i^G, t) \cdot \mathbf{p}_i - (\mathbf{p}_i^T \cdot \nabla\mathbf{v}(\mathbf{x}_i^G, t) \cdot \mathbf{p}_i) \cdot \mathbf{p}_i \quad (13) $$

Again, the simplest integration consists of writing:

$$ \mathbf{p}_i(t + \Delta t) = \mathbf{p}_i(t) + \left( \nabla \mathbf{v}(\mathbf{x}_i^G, t) \cdot \mathbf{p}_i(t) - (\mathbf{p}_i^T(t) \cdot \nabla \mathbf{v}(\mathbf{x}_i^G, t) \cdot \mathbf{p}_i(t)) \cdot \mathbf{p}_i(t) \right) \cdot \Delta t \quad (14) $$

- Compute the stress at position $\mathbf{x}$ and time $t$ by considering all the fibers inside a control volume $\mathcal{V}(\mathbf{x})$ centered at $\mathbf{x}$ by applying

$$ \boldsymbol{\tau}(\mathbf{x}, t) = 2 \cdot \eta \cdot \mathbf{D}(\mathbf{x}, t) + \sum_{i \in \mathcal{V}(\mathbf{x})} \xi \cdot L^2 \cdot ((\nabla \mathbf{v}(\mathbf{x}, t)) : (\mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i)) \quad (15) $$

In the case of coupled models this stress will serve to update the velocity field, however, as previously indicated, in the present paper we will consider that the velocity field is not perturbed by the presence and orientation of the rods.
## 2.2 Mesoscopic description

Because the high number of rods involved in the suspension the description that we just described, and despite its conceptual simplicity, fails to address the situations usually encountered in practice. For this reason coarser descriptions are preferred. The first plausible coarser description applies a zoom-out, in which the rods individuality is loosed in favor of averaged descriptions.

In the case of rods, one could described the microstructure at a certain point $\mathbf{x}$ and time $t$ from the orientation distribution function $\psi(\mathbf{x}, t, \mathbf{p})$ given the fraction of rods that at position $\mathbf{x}$ and time $t$ are oriented in the direction $\mathbf{p}$. Obviously, function $\psi$ verifies the normality condition:

$$ \int_{S(0,1)} \psi(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p} = 1, \quad \forall \mathbf{x}, \forall t \qquad (16) $$

where $S(0, 1)$ is the surface of the unit ball (circumference of unit radius in the 2D case and spherical surface of unit radius in the 3D case) that defines all the possible rod orientation in 2D and 3D respectively.

Thus, we substituted the fine microscopic description, needing for the specification of each rod orientation, by an averaged scalar function. However, in order to use it, one needs to derive the equation governing the evolution of such an orientation distribution function.

The balance ensuring the rods conservation implies:

$$ \frac{\partial \psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\dot{\mathbf{x}} \cdot \psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \cdot \psi) = 0 \qquad (17) $$

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

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \qquad (18) $$

Eq. (17), known as Fokker-Planck equation, is a well balanced compromise between the macroscopic scale that defines the overall process, where the space and time coordinates are defined, and a finer microscopic information describing the rods flow induced orientation mechanism described by the Jeffery's equation (18).

The price to be paid is the increase of the model dimensionality, because the orientation distribution is defined in a high-dimensional domain consisting of 5 dimensions in the general 3D case, i.e. $\psi : (\mathbf{x}, t, \mathbf{p}) \to \mathcal{R}^+$ where $\mathbf{x} \in \Omega \subset \mathcal{R}^3, t \in \mathcal{I} \subset \mathcal{R}^+, \mathbf{p} \in S(0, 1)$.

The extra-stress due to the presence of such a population of rods is determined by averaging the discrete virial stress defined at the microscopic scale

$$ \begin{aligned} \boldsymbol{\tau} &= \boldsymbol{\tau}^f + \boldsymbol{\tau}^r = \\ &= 2 \cdot \eta \cdot \mathbf{D} + 2 \cdot \eta \cdot N_p \cdot \int_{S(0,1)} (\nabla \mathbf{v}) : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}) \cdot \psi(\mathbf{x}, t, \mathbf{p}) \cdot d\mathbf{p} = \\ &= 2 \cdot \eta \cdot \mathbf{D} + 2 \cdot \eta \cdot N_p \cdot (\nabla \mathbf{v} : \mathbf{a}^{(4)}) \end{aligned} \qquad (19) $$

where $\mathbf{a}^{(4)}$ represents the fourth order moment of the orientation distribution function and $N_p$ is the so-called particle number that depends on the rods concentration and
shape. Because the symmetry of $\mathbf{a}^{(4)}$ the gradient of velocities in the previous expression can be replaced by the rate of strain tensor $\mathbf{D}$:

$$ \tau = \tau^f + \tau^r = 2 \cdot \eta \cdot \mathbf{D} + 2 \cdot \eta \cdot N_p \cdot (\mathbf{D} : \mathbf{a}^{(4)}) \quad (20) $$

In summary, the mesoscopic decoupled modeling, consists of:

- Given the flow velocity field $\mathbf{v}(\mathbf{x}, t)$ and the initial orientation distribution $\psi(\mathbf{x}, t = 0, \mathbf{p})$

- Solve the Fokker-Planck equation that governs the rods orientation distribution:

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

with

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (22) $$

- Compute the stress tensor:

$$ \tau = \tau^f + \tau^r = 2 \cdot \eta \cdot \mathbf{D} + 2 \cdot \eta \cdot N_p \cdot (\mathbf{D} : \mathbf{a}^{(4)}) \quad (23) $$

with

$$ \mathbf{a}^{(4)} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \cdot \psi(\mathbf{x}, t, \mathbf{p}) \cdot d\mathbf{p} \quad (24) $$

## 2.3 Macroscopic description

Fokker-Planck based descriptions were rarely considered precisely because the curse of dimensionality that the introduction of the conformation coordinates (the rod orientation in the case here considered) 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 this reason, mesoscopic models were coarsened one step further to derive macroscopic models defined in standard physical domains, involving space and time.

In this section we are illustrating the transition from the mesoscopic to the macroscopic scale. At the macroscopic scale the orientation distribution function is substituted by its moments for describing the microstructure.

The first moment $\mathbf{a}^{(1)}$ results:

$$ \mathbf{a}^{(1)} = \int_{S(0,1)} \mathbf{p} \cdot \psi \, d\mathbf{p} \quad (25) $$

that represents the mean local orientation. Its time evolution can be obtained by taking the time derivative of Eq. (25) and using the expression of $\dot{\mathbf{p}}$ previously deduced:

$$ \begin{aligned} \dot{\mathbf{a}}^{(1)} &= \int_{S(0,1)} \dot{\mathbf{p}} \cdot \psi \, d\mathbf{p} = \int_{S(0,1)} (\nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p}) \cdot \psi \, d\mathbf{p} = \\ &= \nabla \mathbf{v} \cdot \mathbf{a}^{(1)} - \mathbf{a}^{(3)} : \nabla \mathbf{v} \end{aligned} \quad (26) $$
where " : " represents the tensor product twice contracted.

We can notice that the evolution of the first moment involves the third one and so on, needing for a closure relation.

Now, considering the second moment $\mathbf{a}^{(2)}$:

$$ \mathbf{a}^{(2)} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \psi \, d\mathbf{p} \qquad (27) $$

its time derivative reads:

$$ \begin{aligned} \dot{\mathbf{a}}^{(2)} &= \int_{S(0,1)} (\dot{\mathbf{p}} \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}) \cdot \psi \, d\mathbf{p} = \\ &= \int_{S(0,1)} (\nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p}) \otimes \mathbf{p} \cdot \psi \, d\mathbf{p} + \\ &\quad + \int_{S(0,1)} \mathbf{p} \otimes (\nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p}) \cdot \psi \, d\mathbf{p} = \\ &= \nabla \mathbf{v} \cdot \mathbf{a}^{(2)} + \mathbf{a}^{(2)} \cdot (\nabla \mathbf{v})^T - 2 \cdot \mathbf{a}^{(4)} : \nabla \mathbf{v} \end{aligned} \qquad (28) $$

where

$$ \mathbf{a}^{(4)} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \cdot \psi \, d\mathbf{p} \qquad (29) $$

Again a closure relation is needed in order to express the fourth order moment as a function of the lower order moments.

Different closure relations have been introduced and widely used [3] [40]. One of the most used is the quadratic one that reads:

$$ \mathbf{a}^{(4)} = \mathbf{a}^{(2)} \otimes \mathbf{a}^{(2)} \qquad (30) $$

that results exact when the microstructure is locally perfectly aligned in a direction $\hat{\mathbf{p}}(\mathbf{x})$, because in that case the distribution function writes $\psi(\mathbf{x}, t, \mathbf{p}) = \delta(\mathbf{p} - \hat{\mathbf{p}})$ ($\delta(\bullet)$ being the Dirac's delta distribution) and consequently

$$ \hat{\mathbf{a}}^{(4)} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \cdot \delta(\mathbf{p} - \hat{\mathbf{p}}) \, d\mathbf{p} = \hat{\mathbf{p}} \otimes \hat{\mathbf{p}} \otimes \hat{\mathbf{p}} \otimes \hat{\mathbf{p}} = \hat{\mathbf{a}}^{(2)} \otimes \hat{\mathbf{a}}^{(2)} \qquad (31) $$

Moreover, as the model until now described does not involve Brownian effects it can be proved that an exact closure relation exists, known as natural closure relation [29] [52].

In summary, the decoupled macroscopic description consists of:

- Given the flow velocity field $\mathbf{v}(\mathbf{x}, t)$ and the initial second order moment of the orientation distribution $\mathbf{a}^{(2)}(\mathbf{x}, t = 0)$

- Solve the equation governing the evolution of $\mathbf{a}^{(2)}$:

$$ \dot{\mathbf{a}}^{(2)} = \nabla \mathbf{v} \cdot \mathbf{a}^{(2)} + \mathbf{a}^{(2)} \cdot (\nabla \mathbf{v})^T - 2 \cdot \mathbf{a}^{(4)} : \nabla \mathbf{v} \qquad (32) $$

by using an exact or an approximated closure relation.
– Compute the stress tensor:

$$ \tau = \tau^f + \tau^r = 2 \cdot \eta \cdot \mathbf{D} + 2 \cdot \eta \cdot N_p \cdot \left( \mathbf{D} : \mathbf{a}^{(4)} \right) \quad (33) $$

using again an appropriate closure relation, that in the present case, because there is not Brownian effects, an exact closure relation exists.

### 2.3.1 On the objectivity of $\dot{\mathbf{a}}^{(2)}$

One must be careful with the time derivatives of tensors. In our case, when considering the time derivative of the second order moment, and introducing the decomposition of $\nabla \mathbf{v}$ in its symmetric $\mathbf{D}$ and skew-symmetric part $\mathbf{\Omega}$ it results

$$ \dot{\mathbf{a}}^{(2)} = \mathbf{\Omega} \cdot \mathbf{a}^{(2)} - \mathbf{a}^{(2)} \cdot \mathbf{\Omega} + \mathbf{D} \cdot \mathbf{a}^{(2)} + \mathbf{a}^{(2)} \cdot \mathbf{D} - 2 \cdot \mathbf{a}^{(4)} : \mathbf{D} \quad (34) $$

that can be rewritten as

$$ \frac{\delta \mathbf{a}^{(2)}}{\delta t} = \dot{\mathbf{a}}^{(2)} - \mathbf{\Omega} \cdot \mathbf{a}^{(2)} + \mathbf{a}^{(2)} \cdot \mathbf{\Omega} = \mathbf{D} \cdot \mathbf{a}^{(2)} + \mathbf{a}^{(2)} \cdot \mathbf{D} - 2 \cdot \mathbf{a}^{(4)} : \mathbf{D} \quad (35) $$

where we can identify in the left hand member the Oldroyd's derivative $\frac{\delta \mathbf{a}^{(2)}}{\delta t}$ and in the right hand member only dependences on the rate of strain tensor. Form that we can conclude on the objectivity of the equation governing the time evolution of the second order moment $\mathbf{a}^{(2)}$.

In what follows and for the sake of clarity we will denote by $\mathbf{a}$ the second order moment $\mathbf{a}^{(2)}$ and by $\mathbf{A}$ the fourth order one $\mathbf{a}^{(4)}$.

## 3 Dilute suspensions of Brownian rods

Until now Brownian effects were neglected. The Brownian effects are due to the fluid molecules bombardment acting on the beads. These effects were widely analyzed in [32] [33] [25] when focusing in a microscopic description.

### 3.1 Microscopic description

In this case the road beads are subjected to the hydrodynamical forces and the ones coming from such bombardment. The first one was introduced previously

$$ \mathbf{F}^H = \xi \cdot (\nabla \mathbf{v} \cdot \mathbf{p} \cdot L - \dot{\mathbf{p}} \cdot L) \quad (36) $$

where the superscript "H" refers to its hydrodynamic nature. Now the Brownian force $\mathbf{F}_I^B$ is assumed applying during a short time interval $\delta t$ following a certain statistical distribution concerning its magnitude and its orientation. The first one is assumed described by a gaussian distribution of zero mean and a certain standard deviation and the one related to the orientation by a uniform distribution on the unit circle (2D) or in the unit sphere (3D).

Brownian forces applying in the rod direction are assumed equilibrated. However, the components of those forces perpendicular to the rod axis contribute to the rod rotation, and then they affect the rod rotational velocity. In what follows, for the sake
of clarity, we restrict our analysis to the 2D case. Because the rod inertialess, the resultant moment must vanish:

$$ \mathbf{F}^H \cdot \mathbf{t} + \mathbf{F}^B \cdot \mathbf{t} = 0 \qquad (37) $$

being $\mathbf{t}$ the unit vector tangent to the unit cercle. By introducing in this balance the expression of the hydrodynamic force it results

$$ \mathbf{t}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p} - \|\dot{\mathbf{p}}\| = -\frac{\mathbf{F}^B \cdot \mathbf{t}}{\xi \cdot L} \qquad (38) $$

and taking into account that necessarily $\dot{\mathbf{p}} = \|\dot{\mathbf{p}}\| \cdot \mathbf{t}$, it results

$$ \dot{\mathbf{p}} = (\mathbf{t}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{t} + \frac{\mathbf{F}^B \cdot \mathbf{t}}{\xi \cdot L} \cdot \mathbf{t} \qquad (39) $$

that using the vectorial equivalence

$$ (\mathbf{t}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{t} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \qquad (40) $$

results in

$$ \begin{aligned} \dot{\mathbf{p}} &= \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} + \frac{\mathbf{F}^B \cdot \mathbf{t}}{\xi \cdot L} \cdot \mathbf{t} = \\ &= \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} + \frac{\mathbf{F}^B - (\mathbf{F}^B \cdot \mathbf{p}) \cdot \mathbf{p}}{\xi \cdot L} \end{aligned} \qquad (41) $$

where we can notice that the rotation velocity is given by the Jeffery expression $\dot{\mathbf{p}}^J$ complemented with a term accounting for the Bownian effects $\dot{\mathbf{p}}^B$:

$$ \dot{\mathbf{p}}^J = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \qquad (42) $$

and

$$ \dot{\mathbf{p}}^B = \frac{\mathbf{F}^B - (\mathbf{F}^B \cdot \mathbf{p}) \cdot \mathbf{p}}{\xi \cdot L} \qquad (43) $$

from which

$$ \dot{\mathbf{p}} = \dot{\mathbf{p}}^J + \dot{\mathbf{p}}^B \qquad (44) $$

Finally, we must discuss on the effects of such Brownian contribution to the extra-stress tensor.

Consider a rod aligned along the $x$-axis, such that $\mathbf{p}^T = (1, 0)$ and the fluid at rest. This rod is subjected to a continuous bombardment from the solvent molecules (remember that the regime is dilute enough for neglecting the rod-rod interactions). The components of forces aligned with the rod axis do not participate to the rod rotation and by averaging on a time $\Delta t \gg \delta t$ the contributions of the components of those forces along the rod direction vanish. On the contrary, the ones perpendicular to the rod will participate to the stress. To derive the expression of this Brownian contribution we consider that due to a Brownian force the rod rotates a small angle $\delta\theta > 0$, being the rod orientation defined by $\mathbf{p}_{\delta\theta}$. Considering the Brownian force applying at that position [32] [33], i.e. $\|\mathbf{F}^B\| \cdot \mathbf{t}_{\delta\theta}$, it results the contribution to the virial stress given by:
$$ -\|\mathbf{F}^B\| \cdot \mathbf{p}_{\delta\theta} \otimes \mathbf{t}_{\delta\theta} = \|\mathbf{F}^B\| \cdot \begin{pmatrix} \sin(\delta\theta) \cdot \cos(\delta\theta) & -\cos^2(\delta\theta) \\ \sin^2(\delta\theta) & -\sin(\delta\theta) \cdot \cos(\delta\theta) \end{pmatrix} \quad (45) $$

where the negative sign accounts for the fact that the hydrodynamic force applies in the opposite direction of the Brownian force.

We can notice in that expression two facts: (i) the trace is zero, and (ii) the contribution is non-symmetric. However, we can imagine that a little bit later, a Brownian force will apply in the opposite direction, leading to an angle $-\delta\theta$, from which

$$ \|\mathbf{F}^B\| \cdot \mathbf{p}_{-\delta\theta} \otimes \mathbf{t}_{-\delta\theta} = \|\mathbf{F}^B\| \cdot \begin{pmatrix} \sin(\delta\theta) \cdot \cos(\delta\theta) & \cos^2(\delta\theta) \\ -\sin^2(\delta\theta) & -\sin(\delta\theta) \cdot \cos(\delta\theta) \end{pmatrix} \quad (46) $$

and then after averaging in $\Delta t$ it results a Brownian contribution to the extra-stress due to the rod aligned on the x-direction ($\varphi = 0$)

$$ \tau_{\varphi=0}^{B} \approx \begin{pmatrix} \beta & 0 \\ 0 & -\beta \end{pmatrix} = \beta \cdot \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} = \beta \cdot \mathbf{U} \quad (47) $$

that becomes almost symmetric and traceless.

Now, for rods aligned in any other direction $\varphi$, it suffices to apply a rotation of angle $\varphi$ to tensor $\tau_{\varphi=0}^B$:

$$ \tau_{\varphi}^{B} = \beta \cdot \mathbf{R}_{\varphi}^{T} \cdot \mathbf{U} \cdot \mathbf{R}_{\varphi} \quad (48) $$

with

$$ R_{\varphi} = \begin{pmatrix} \cos(\varphi) & \sin(\varphi) \\ -\sin(\varphi) & \cos(\varphi) \end{pmatrix} \quad (49) $$

that finally results:

$$ \tau_{\varphi}^{B} = 2 \cdot \beta \cdot \left( \begin{pmatrix} \cos^2(\varphi) & \sin(\varphi) \cdot \cos(\varphi) \\ \sin(\varphi) \cdot \cos(\varphi) & \sin^2(\varphi) \end{pmatrix} - \begin{pmatrix} \frac{1}{2} & 0 \\ 0 & \frac{1}{2} \end{pmatrix} \right) \quad (50) $$

that can be written in compact form

$$ \tau_{\mathbf{p}_i}^B = 2 \cdot \beta \cdot \left( \mathbf{p}_i \otimes \mathbf{p}_i - \frac{\mathbf{I}}{2} \right) \quad (51) $$

where $\mathbf{I}$ is the identity tensor.

For a population of rods $\mathbf{p}_i$, $i = 1, \dots, N$, the contribution of Brownian effects results finally:

$$ \tau^B = 2 \cdot \beta \cdot \sum_{i=1}^{i=N} \left( \mathbf{p}_i \otimes \mathbf{p}_i - \frac{\mathbf{I}}{2} \right) \quad (52) $$

The microscopic description can be summarized as follows:

- Given the flow velocity field $\mathbf{v}(\mathbf{x}, t)$, the initial position of the centre of gravity of each rod $\mathbf{x}_i^G(t = 0)$ and its orientation $\mathbf{p}_i(t = 0)$ or knowing both them at time $t$

- For each fiber $i = 1, \dots, N$:
- Update the position of the centre of gravity by integrating:

$$ \frac{d\mathbf{x}_i^G}{dt} = \mathbf{v}(\mathbf{x}_i^G, t) \qquad (53) $$

- Update its orientation by integrating

$$ \frac{d\mathbf{p}_i}{dt} = \nabla \mathbf{v}(\mathbf{x}_i^G, t) \cdot \mathbf{p}_i - (\mathbf{p}_i^T \cdot \nabla \mathbf{v}(\mathbf{x}_i^G, t) \cdot \mathbf{p}_i) \cdot \mathbf{p}_i + \frac{\mathbf{F}_i^B - (\mathbf{F}_i^B \cdot \mathbf{p}_i) \cdot \mathbf{p}_i}{\xi \cdot L} \qquad (54) $$

- Compute the stress at position x and time t by considering all the fibers inside a control volume V(x) centered at x by applying

$$ \tau(\mathbf{x}, t) = 2 \cdot \eta \cdot \mathbf{D}(\mathbf{x}, t) + \sum_{i \in \mathcal{V}(\mathbf{x})} \left( \xi \cdot L^2 \cdot \left( \left( \nabla \mathbf{v}(\mathbf{x}, t) \right) : \left( \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \right) \right) + 2 \cdot \beta \cdot \left( \mathbf{p}_i \otimes \mathbf{p}_i - \frac{\mathbf{I}}{2} \right) \right) \qquad (55) $$

In the case of coupled models this stress will serve to update the velocity field, however, as previously indicated, in the present paper we will consider that the velocity field is not perturbed by the presence and orientation of the rods.

## 3.2 Mesoscopic description

At the mesoscopic scale we postulate that Brownian effects try to randomize the rods orientation distribution, mechanism that can be modeled by assuming a diffusion term in the Fokker-Planck equation:

$$ \frac{\partial \psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \cdot \psi) = D_r \frac{\partial^2 \psi}{\partial \mathbf{p}^2} \qquad (56) $$

with the same flow induced orientation term given by the Jeffery's equation:

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \qquad (57) $$

and where $D_r$ is the rotary diffusion.

We can notice that in absence of flow, i.e. $\mathbf{v}(\mathbf{x}, t) = \mathbf{0}$ the Fokker-Planck equation reduces to

$$ \frac{\partial \psi}{\partial t} = D_r \frac{\partial^2 \psi}{\partial \mathbf{p}^2} \qquad (58) $$

that ensures an steady state isotropic orientation distribution: $\psi(\mathbf{x}, t \to \infty, \mathbf{p}) = \frac{1}{2 \cdot \pi}$ in 2D and $\psi(\mathbf{x}, t \to \infty, \mathbf{p}) = \frac{1}{4 \cdot \pi}$ in 3D. The higher is the rotational diffusion the faster the isotropic orientation distribution in reached.

Thus, at the mesoscopic level the introduction of Brownian effects seems quite simple. The question is what is the microscopic and macroscopic counterpart of the just introduced diffusion term?

The Fokker-Planck equation can be rewritten in the form:

$$ \frac{\partial \psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \cdot \psi) - D_r \frac{\partial^2 \psi}{\partial \mathbf{p}^2} = $$
$$ = \frac{\partial \psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \cdot \psi) = 0 \qquad (59) $$

where the effective rotational velocity $\dot{\mathbf{p}}$ is given by

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} - D_r \frac{\frac{\partial \psi}{\partial \mathbf{p}}}{\psi} \qquad (60) $$

which contains the flow induced Jeffery's contribution $\dot{\mathbf{p}}^J$ plus the Brownian one $\dot{\mathbf{p}}^B$, i.e. $\dot{\mathbf{p}} = \dot{\mathbf{p}}^J + \dot{\mathbf{p}}^B$ with

$$ \dot{\mathbf{p}}^J = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \qquad (61) $$

and

$$ \dot{\mathbf{p}}^B = -D_r \frac{\frac{\partial \psi}{\partial \mathbf{p}}}{\psi} \qquad (62) $$

The Fokker-Planck multidimensionality issue was usually circumvented by using stochastic strategies on the Ito's equation related to its Fokker-Planck counterpart [50]. It is important to mention that any Fokker-Planck equation has an Ito's counterpart, however, not all the Ito's stochastic equations has necessarily a Fokker-Planck counterpart.

It is also important to mention that the Brownian rotational velocity just defined has not a physical meaning because as we are considering a dilute enough suspension, rod-rod interactions are excluded and then the expression relating Brownian velocity with the gradient of the orientation distribution is meaningfulness.

At the mesoscopic scale the Brownian contribution to the extra-stress tensor results from the generalization of Eq. (52):

$$ \boldsymbol{\tau}^B = 2 \cdot \gamma \cdot \int_{S(0,1)} \left( \mathbf{p} \otimes \mathbf{p} - \frac{\mathbf{I}}{2} \right) \cdot \psi(\mathbf{p}) \, d\mathbf{p} = 2 \cdot \gamma \left( \mathbf{a} - \frac{\mathbf{I}}{2} \right) \qquad (63) $$

In summary, the mesoscopic decoupled description consists of:

- Given the flow velocity field $\mathbf{v}(\mathbf{x}, t)$ and the initial orientation distribution $\psi(\mathbf{x}, t = 0, \mathbf{p})$

- Solve the Fokker-Planck equation that governs the rods orientation distribution:

$$ \frac{\partial \psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \cdot \psi) = D_r \cdot \frac{\partial^2 \psi}{\partial \mathbf{p}^2} \qquad (64) $$

with

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \qquad (65) $$

- Compute the stress tensor:

$$ \boldsymbol{\tau} = \boldsymbol{\tau}^f + \boldsymbol{\tau}^r = 2 \cdot \eta \cdot \mathbf{D} + 2 \cdot \eta \cdot N_p \cdot (\mathbf{D} : \mathbf{A}) + 2 \cdot \gamma \cdot \left( \mathbf{a} - \frac{\mathbf{I}}{2} \right) \qquad (66) $$
### 3.3 Macroscopic description

When moving towards the macroscopic scale the Brownian contribution to the extra-stress is defined by Eq. (63), however at the macroscopic scale the microstructure is defined by the different moments of the orientation distribution. In what follows we are deriving the contribution of Brownian effects on the equation governing the evolution of the second order moment.

We start from the second order moment definition

$$ \mathbf{a} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \psi \, d\mathbf{p} \qquad (67) $$

whose time derivative involves now the effective rotational velocity $\dot{\mathbf{p}}$

$$ \dot{\mathbf{a}} = \int_{S(0,1)} (\dot{\mathbf{p}} \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}) \cdot \psi \, d\mathbf{p} \qquad (68) $$

As proved in the previous sections, the effective rotational velocity contains the flow induced contribution given by the Jeffery's expression $\dot{\mathbf{p}}^J$ and the one induced by the Brownian effects $\dot{\mathbf{p}}^B$, given by Eqs. (61) and (62) respectively. With this decomposition Eq. (68) can be written as

$$ \begin{aligned} \dot{\mathbf{a}} = \int_{S(0,1)} & \left( (\dot{\mathbf{p}}^J + \dot{\mathbf{p}}^B) \otimes \mathbf{p} + \mathbf{p} \otimes (\dot{\mathbf{p}}^J + \dot{\mathbf{p}}^B) \right) \cdot \psi \, d\mathbf{p} = \\ & = \int_{S(0,1)} \left( \dot{\mathbf{p}}^J \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}^J \right) \cdot \psi \, d\mathbf{p} + \\ & + \int_{S(0,1)} \left( \dot{\mathbf{p}}^B \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}^B \right) \cdot \psi \, d\mathbf{p} = \dot{\mathbf{a}}^J + \dot{\mathbf{a}}^B \end{aligned} \qquad (69) $$

where the flow induced microstructure evolution $\dot{\mathbf{a}}^J$ is given by Eq. (32)

$$ \dot{\mathbf{a}}^J = \nabla \mathbf{v} \cdot \mathbf{a} + \mathbf{a} \cdot (\nabla \mathbf{v})^T - 2 \cdot \mathbf{A} : \mathbf{D} \qquad (70) $$

We are now calculating the expression of the remaining contribution $\dot{\mathbf{a}}^B$:

$$ \dot{\mathbf{a}}^B = \int_{S(0,1)} \left( \dot{\mathbf{p}}^B \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}^B \right) \cdot \psi \, d\mathbf{p} \qquad (71) $$

with $\dot{\mathbf{p}}^B$ given by

$$ \dot{\mathbf{p}}^B = -D_r \frac{\frac{\partial \psi}{\partial \mathbf{p}}}{\psi} \qquad (72) $$

For the sake of clarity we consider again the 2D case that allows writing

$$ \dot{\mathbf{p}}^B = -D_r \frac{\frac{\partial \psi}{\partial \theta}}{\psi} \cdot \mathbf{t} \qquad (73) $$

with $\mathbf{t}$ being the unit tangent vector to the unit circle. In this case Eq. (71) reduces to

$$ \dot{\mathbf{a}}^B = -D_r \cdot \int_{S(0,1)} (\mathbf{t} \otimes \mathbf{p} + \mathbf{p} \otimes \mathbf{t}) \cdot \frac{\partial \psi}{\partial \theta} \, d\theta \qquad (74) $$
Now, integrating by parts Eq. (74) and taking into account

$$ \frac{d\mathbf{p}}{d\theta} = \mathbf{t} \qquad (75) $$

and

$$ \frac{dt}{d\theta} = -\mathbf{p} \qquad (76) $$

it results:

$$ \dot{\mathbf{a}}^B = -2 \cdot D_r \cdot \int_{S(0,1)} (\mathbf{p} \otimes \mathbf{p} - \mathbf{t} \otimes \mathbf{t}) \cdot \psi(\theta) \, d\theta \qquad (77) $$

It is easy to prove that

$$ \mathbf{t} \otimes \mathbf{t} + \mathbf{p} \otimes \mathbf{p} = \mathbf{I} \quad \rightarrow \quad \mathbf{t} \otimes \mathbf{t} = \mathbf{I} - \mathbf{p} \otimes \mathbf{p} \qquad (78) $$

that allows writing Eq. (77) in the form

$$ \dot{\mathbf{a}}^B = -2 \cdot D_r \cdot \int_{S(0,1)} (2 \cdot \mathbf{p} \otimes \mathbf{p} - \mathbf{I}) \cdot \psi(\theta) \, d\theta = -2 \cdot D_r \cdot (2\mathbf{a} - \mathbf{I}) \qquad (79) $$

or

$$ \dot{\mathbf{a}}^B = -4 \cdot D_r \cdot \left( \mathbf{a} - \frac{\mathbf{I}}{2} \right) \qquad (80) $$

We can notice that in absence of flow, $\dot{\mathbf{a}}^J = \mathbf{0}$, and then $\dot{\mathbf{a}} = \dot{\mathbf{a}}^B$

$$ \dot{\mathbf{a}} = -4 \cdot D_r \cdot \left( \mathbf{a} - \frac{\mathbf{I}}{2} \right) \qquad (81) $$

ensuring an isotropic steady state, i.e. $\mathbf{a}(t \to \infty) = \frac{\mathbf{I}}{2}$ in the 2D case.

In summary, the decoupled macroscopic description consists of:

- Given the flow velocity field $\mathbf{v}(\mathbf{x}, t)$ and the initial second order moment of the orientation distribution $\mathbf{a}(\mathbf{x}, t = 0)$

- Solve the equation governing the evolution of $\mathbf{a}$:

$$ \dot{\mathbf{a}} = \nabla \mathbf{v} \cdot \mathbf{a} + \mathbf{a} \cdot (\nabla \mathbf{v})^T - 2 \cdot \mathbf{A} : \mathbf{D} - 4 \cdot D_r \cdot \left( \mathbf{a} - \frac{\mathbf{I}}{2} \right) \qquad (82) $$

by using an exact or an approximated closure relation.

- Compute the stress tensor:

$$ \boldsymbol{\tau} = \boldsymbol{\tau}^f + \boldsymbol{\tau}^r = 2 \cdot \eta \cdot \mathbf{D} + 2 \cdot \eta \cdot N_p \cdot (\mathbf{D} : \mathbf{A}) + 2 \cdot \gamma \cdot \left( \mathbf{a} - \frac{\mathbf{I}}{2} \right) \qquad (83) $$

using again an appropriate closure relation, that in the present case due to the presence of a diffusion term will be no more exact.
## 4 Semi-concentrated regime

Semi-dilute and semi-concentrated regimes have been widely addressed, most of time by using phenomenological approaches. The most common approach consists of considering that rod-rod interactions tends to randomize the orientation distribution. Thus, a second diffusion coefficient is introduced for accounting for rods interactions, however, in the present case that diffusion coefficient should scale with the flow intensity in order to ensure that in absence of flow the microstructure does not evolve artificially because such a diffusion term. In general the interaction diffusion coefficient $D_I$ is assumed in the general form

$$D_I = C_I \cdot f(D^{eq}) \qquad (84)$$

where $D^{eq}$ is related to the second invariant of the rate of strain tensor, i.e. $D^{eq} = \sqrt{\overline{\mathbf{D}} : \overline{\mathbf{D}}}$. [31] considered the simplest dependence $f(D^{eq}) = D^{eq}$.

Obviously, there are finer approaches based on the direct simulation where the rod-rod interactions are taken into account explicitly [28]. A nice mesoscopic based macroscopic approach was proposed in [30].

In what follows we are considering a different approach based on the kinematics and dynamics of entangled rods generating aggregates or clusters, allowing to derive an, in our knowledge original two-scales kinetic theory approach.

## 5 Suspensions involving rigid clusters

As just discussed, when the concentration is high enough, a sort of clusters composed of entangled rods are observed. When the suspension flows, those clusters seem animated of almost rigid motions.

Clusters can be sketched in two dimensions as entangled aggregates of rods (see Fig. 3). As in the previous sections rods can be idealized by a rigid segment joining the two opposite beads. Forces are assumed applying on the beads. Microscopic simulations become unaffordable because the computational resources needed to address scenarios of real interest. On the other extremum, macroscopic modeling is without any doubt too phenomenological for describing the main flow and microstructure features. Kinetic theory approaches could be again an appealing compromise between accuracy and computational efficiency.

To define the cluster representation that we will consider in the analysis that follows we extract the beads located on the periphery of the cluster as illustrated in Fig. 4. Then we consider the star configuration depicted in Fig. 5 that will be the cluster idealization on which we are performing the analysis of its kinematics and dynamics.

Obviously, this process constitutes a tremendous simplification but it constitutes the simplest model that could be enriched without major difficulties to cover more complex and rich systems.

### 5.1 Microscopic description of rigid clusters composed of rods

First we are considering a rigid cluster consisting of N rods, oriented in the directions $\mathbf{p}_i$ as sketched in Fig. 6. Brownian effects can be neglected and then only flow induced
<figure><img src="image_4.png" /><figcaption>Fig. 3 Entangled rods constituting a rigid cluster</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Extracting the active beads located on the cluster periphery</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Defining the rigid star representation</figcaption></figure>

hydrodynamic forces must be considered. A force $\mathbf{F}_i$ applies on each bead $L \cdot \mathbf{p}_i$, assumed again given by:

$$ \mathbf{F}_i = \xi \cdot (\nabla \mathbf{v} \cdot \mathbf{p}_i \cdot L - \dot{\mathbf{p}}_i \cdot L) \quad (85) $$

Now, instead of enforcing a null resulting moment for each rod as was the case when addressing dilute solutions, we are enforcing a null resulting moment for the whole cluster.

The moment associated with rod $i$ is given by
<figure><img src="image_7.png" /><figcaption>Fig. 6 Hydrodynamic forces applying on a rigid cluster composed of rods</figcaption></figure>

$$ M_i = 2 \cdot L \cdot \mathbf{p}_i \times \mathbf{F}_i \qquad (86) $$

where the symbol " × " denotes the cross product.

The forces are by construction self-equilibrated for the whole cluster, thus we only need to enforce the nullity of the resulting moment:

$$ \sum_{i=1}^{i=N} M_i = 0 \qquad (87) $$

that taking into account Eqs. (86) and (85) results in

$$ \sum_{i=1}^{i=N} \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i) = \sum_{i=1}^{i=N} \mathbf{p}_i \times \dot{\mathbf{p}}_i \qquad (88) $$

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

$$ \dot{\mathbf{p}}_i = \omega \times \mathbf{p}_i \qquad (89) $$

then:

- In the 2D case, and because $\omega$ and $\mathbf{p}_i$ are orthogonal, $\mathbf{p}_i \times \dot{\mathbf{p}}_i = \omega \forall i$, implying:

$$ \sum_{i=1}^{i=N} \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i) = \sum_{i=1}^{i=N} \mathbf{p}_i \times \dot{\mathbf{p}}_i = N \cdot \omega \qquad (90) $$

Thus, the kinematics of the rigid 2D cluster can be defined from:

$$ \omega = \frac{\sum_{i=1}^{i=N} \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i)}{N} \qquad (91) $$

being the angular velocity of each rod j:

$$ \dot{\mathbf{p}}_j = \omega \times \mathbf{p}_j = \frac{\left(\sum_{i=1}^{i=N} \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i)\right) \times \mathbf{p}_j}{N} \qquad (92) $$
Remark: In the case of a single rod $N = 1$ whose orientation is defined by **p**, the previous expression reduces to

$$ \dot{\mathbf{p}} = \omega \times \mathbf{p} = (\mathbf{p} \times (\nabla \mathbf{v} \cdot \mathbf{p})) \times \mathbf{p} \quad (93) $$

and using the vector triple product formula $\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = \mathbf{b} \cdot (\mathbf{a} \cdot \mathbf{c}) - \mathbf{c} \cdot (\mathbf{a} \cdot \mathbf{b})$ it results

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (94) $$

that corresponds with the Jeffery's equation for rods. □

– In the general 3D case, $\omega$ and $\mathbf{p}_i$ are not necessarily orthogonal and then

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

and using again the vector triple product formula it results:

$$ \mathbf{p}_i \times \dot{\mathbf{p}}_i = \mathbf{p}_i \times (\omega \times \mathbf{p}_i) = \omega - (\omega \cdot \mathbf{p}_i) \cdot \mathbf{p}_i \quad (96) $$

Thus, we obtain in this case from (88)

$$ \sum_{i=1}^{i=N} \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i) = N \cdot \omega - \sum_{i=1}^{i=N} (\omega \cdot \mathbf{p}_i) \cdot \mathbf{p}_i = N \cdot \omega - \sum_{i=1}^{i=N} (\mathbf{p}_i \otimes \mathbf{p}_i) \cdot \omega \quad (97) $$

that can be rewritten in a more compact form as:

$$ \omega = \left( N \cdot \mathbf{I} - \sum_{i=1}^{i=N} (\mathbf{p}_i \otimes \mathbf{p}_i) \right)^{-1} \cdot \left( \sum_{i=1}^{i=N} \mathbf{p}_i \times (\nabla \mathbf{v} \cdot \mathbf{p}_i) \right) \quad (98) $$

## 5.2 Jeffery revisited

Until now we addressed the Jeffery's solution for infinite aspect ratio ellipsoids, i.e., slender bodies or rods. The 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 [38]:

$$ \dot{\mathbf{p}} = \mathbf{\Omega} \cdot \mathbf{p} + r \cdot \mathbf{D} \cdot \mathbf{p} - r (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (99) $$

where $r$ depends on the ellipsoid aspect ratio $\lambda$, the former given by the length-diameter ratio:

$$ r = \frac{\lambda^2 - 1}{\lambda^2 + 1} \quad (100) $$

It can be noticed that by considering $\lambda \to \infty$, $r \approx 1$ and then Eq. (99) reduces to

$$ \dot{\mathbf{p}} = (\mathbf{\Omega} + \mathbf{D}) \cdot \mathbf{p} - (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (101) $$

or

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (102) $$
<figure><img src="image_8.png" /><figcaption>Fig. 7 Two-rods rigid cluster representing an ellipsoidal 2D particle</figcaption></figure>

by considering the equality $\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p} = \mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}$. This expression corresponds to the expression previously derived by considering a simple momentum balance of a two-beads rod.

We are proving that in order to represents non-infinite aspect ratio 2D particles it suffices to consider a rigid system composed of two rods, aligning perpendicularly the one with respect to the other, and having lengths $2 \cdot L_1$ and $2 \cdot L_2$, $L_1 \ge L_2$, as depicted in Fig. 7, such that the parameter $r$ in the general Jeffery's equation (99) results from $\lambda = \frac{L_1}{L_2}$

$$r = \frac{L_1^2 - L_2^2}{L_1^2 + L_2^2} \qquad (103)$$

In the present configuration the forces applying in bead $L_1 \cdot \mathbf{p}_1$ and $L_2 \cdot \mathbf{p}_2$ read

$$\mathbf{F}_1 = \xi \cdot L_1 \cdot (\nabla \mathbf{v} \cdot \mathbf{p}_1 - \dot{\mathbf{p}}_1) \qquad (104)$$

and

$$\mathbf{F}_2 = \xi \cdot L_2 \cdot (\nabla \mathbf{v} \cdot \mathbf{p}_2 - \dot{\mathbf{p}}_2) \qquad (105)$$

respectively, with

$$\left\{ \begin{array}{l} \dot{\mathbf{p}}_1 = \omega \times \mathbf{p}_1 \\ \dot{\mathbf{p}}_2 = \omega \times \mathbf{p}_2 \end{array} \right. \qquad (106)$$

The moment balance implies:

$$L_1^2 \cdot \mathbf{p}_1 \times (\nabla \mathbf{v} \cdot \mathbf{p}_1 - \dot{\mathbf{p}}_1) + L_2^2 \cdot \mathbf{p}_2 \times (\nabla \mathbf{v} \cdot \mathbf{p}_2 - \dot{\mathbf{p}}_2) = \mathbf{0} \qquad (107)$$

Introducing (106) into (107) 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$, it results

$$\omega = \frac{L_1^2}{L_1^2 + L_2^2} \cdot (\mathbf{p}_1 \times (\nabla \mathbf{v} \cdot \mathbf{p}_1)) + \frac{L_2^2}{L_1^2 + L_2^2} \cdot (\mathbf{p}_2 \times (\nabla \mathbf{v} \cdot \mathbf{p}_2)) \qquad (108)$$
Now, the ellipsoid velocity of rotation will be given by:

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

It is easy to verify that as expected $\dot{\mathbf{p}}_1 \cdot \mathbf{p}_1 = 0$, i.e. the variation of $\mathbf{p}_1$ is directed along its perpendicular direction that 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. (109) reads:

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

We are developing the last term in Eq. (110) in order to obtain an equation that only contains $\mathbf{p}_1$ in order to compare with the Jeffery's expression (99).

First we apply the decomposition $\nabla \mathbf{v} = \mathbf{D} + \mathbf{\Omega}$ from with the last term in Eq. (110) writes:

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

We proceed from the first term in the right hand member of Eq. (111), using the fact that $\mathbf{D}$ is symmetric

$$ 
\begin{aligned}
(\mathbf{p}_1^T \cdot \mathbf{D} \cdot \mathbf{p}_2) \cdot \mathbf{p}_2 &= \\
&= \mathbf{p}_2 \cdot (\mathbf{p}_1^T \cdot \mathbf{D} \cdot \mathbf{p}_2) = \mathbf{p}_2 \cdot (\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 (112) $$

Being $\mathbf{p}_1$ and $\mathbf{p}_2$ mutually perpendicular, it results

$$ (\mathbf{p}_2 \otimes \mathbf{p}_2) + (\mathbf{p}_1 \otimes \mathbf{p}_1) = \mathbf{I} \quad (113) $$

and then Eq. (112) can be written as

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

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

$$ 
\begin{aligned}
(\mathbf{p}_1^T \cdot \mathbf{\Omega} \cdot \mathbf{p}_2) \cdot \mathbf{p}_2 &= \\
&= \mathbf{p}_2 \cdot (\mathbf{p}_1^T \cdot \mathbf{\Omega} \cdot \mathbf{p}_2) = -\mathbf{p}_2 \cdot (\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 (115) $$

that finally results in

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

Now, coming back to Eq. (110), 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} \cdot \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) \cdot \mathbf{p}_1 \quad (117) $$

that corresponds exactly with the Jeffery's expression for an ellipsoid of aspect ratio $\frac{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 (118) $$

### 5.3 Mesoscopic description of rigid clusters composed of rods

When the cluster contains many rods, an alternative description consists of defining its orientation distribution $\psi(\mathbf{p})$.

In this case all the sums in the previous section must be substituted by the corresponding integrals weighted with the distribution function $\psi$.

Thus in the 2D case Eq. (91) leads to:

$$ \omega = \int_{S(0,1)} \mathbf{p} \times (\nabla \mathbf{v} \cdot \mathbf{p}) \cdot \psi(\mathbf{p}) \, d\mathbf{p} \quad (119) $$

being the angular velocity of any rod $\mathbf{q}$:

$$ \dot{\mathbf{q}} = \omega \times \mathbf{q} = \int_{S(0,1)} (\mathbf{p} \times (\nabla \mathbf{v} \cdot \mathbf{p})) \times \mathbf{q} \cdot \psi(\mathbf{p}) \, d\mathbf{p} \quad (120) $$

Cross products can be written in a more compact form by using the permutation of Levi-Civita symbol $\epsilon$ that states $(\mathbf{a} \times \mathbf{b})_i = \epsilon_{ijk} \cdot \mathbf{a}_j \cdot \mathbf{b}_k = \epsilon : (\mathbf{a} \otimes \mathbf{b})$.

Because $(\mathbf{p} \times (\nabla \mathbf{v} \cdot \mathbf{p}))_i = \epsilon_{ijk} \cdot \mathbf{p}_j \cdot (\nabla \mathbf{v})_{kl} \cdot \mathbf{p}_l$, it is easy to prove that Eq. (119) writes:

$$ \omega = \epsilon : (\mathbf{a} \cdot (\nabla \mathbf{v})^T) \quad (121) $$

On the other hand we can define the rotation tensor $\mathbf{W}$ such that

$$ \omega \times \mathbf{p} = \mathbf{W} \cdot \mathbf{p} \quad (122) $$

that allows writing

$$ \dot{\mathbf{p}} = \mathbf{W} \cdot \mathbf{p} \quad (123) $$

As in the 2D case $\mathbf{p}^T = (\mathbf{p}_1, \mathbf{p}_2, 0)$ and $\omega^T = (0, 0, \omega_3)$, the rotation tensor writes:

$$ \mathbf{W} = \begin{pmatrix} 0 & -\omega_3 \\ \omega_3 & 0 \end{pmatrix} \quad (124) $$

Thus, the mesoscopic description uses:

$$ \omega = \epsilon : (\mathbf{a} \cdot (\nabla \mathbf{v})^T) \quad (125) $$

and

$$ \dot{\mathbf{p}} = \omega \times \mathbf{p} \equiv \mathbf{W} \cdot \mathbf{p} \quad (126) $$
with $\mathbf{W} = \mathbf{W}(\omega)$

Thus, a rigid cluster rotates with a velocity that only depends on the second moment of its orientation distribution function $\mathbf{a}$. Any cluster having the same $\mathbf{a}$ will have the same rotational velocity. In this case we should derive an expression given the evolution of $\mathbf{a}$. For this purpose, we start from the definition of $\mathbf{a}$

$$ \mathbf{a} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \psi(\mathbf{p}) \, d\mathbf{p} \qquad (127) $$

and apply the time derivative

$$ \dot{\mathbf{a}} = \int_{S(0,1)} (\dot{\mathbf{p}} \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}) \cdot \psi(\mathbf{p}) \, d\mathbf{p} \qquad (128) $$

Now, taking into account Eq. (126) it results

$$ \begin{aligned} \dot{\mathbf{a}} &= \int_{S(0,1)} ((\mathbf{W} \cdot \mathbf{p}) \otimes \mathbf{p} + \mathbf{p} \otimes (\mathbf{W} \cdot \mathbf{p})) \cdot \psi(\mathbf{p}) \, d\mathbf{p} = \\ &= \mathbf{W} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{W}^T \end{aligned} \qquad (129) $$

### 5.3.1 On the objectivity of $\dot{\mathbf{a}}$

Before going forward we should check the objectivity of expression (129). For this purpose we develop the expression of $\omega_3$ involved in $\mathbf{W}$, obtaining:

$$ \omega_3 = f(\mathbf{D}, \mathbf{a}) - \mathbf{\Omega}_{12} \qquad (130) $$

with

$$ f(\mathbf{D}, \mathbf{a}) = -\mathbf{D}_{12} - 2 \cdot \mathbf{D}_{11} \cdot \mathbf{a}_{12} + 2 \cdot \mathbf{D}_{12} \cdot \mathbf{a}_{11} \qquad (131) $$

Thus $\mathbf{W}$ can be decomposed according to

$$ \mathbf{W} = \begin{pmatrix} 0 & -f \\ f & 0 \end{pmatrix} + \begin{pmatrix} 0 & \mathbf{\Omega}_{12} \\ -\mathbf{\Omega}_{12} & 0 \end{pmatrix} = \mathbf{F} + \mathbf{\Omega} \qquad (132) $$

Thus, Eq. (129) reads

$$ \dot{\mathbf{a}} = \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F} \qquad (133) $$

or

$$ \dot{\mathbf{a}} - \mathbf{\Omega} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{\Omega} = \mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F} \qquad (134) $$

where the left hand member corresponds to the Jauman's objective derivative $\frac{\delta_0 \mathbf{a}}{\delta t}$

$$ \frac{\delta_0 \mathbf{a}}{\delta t} \equiv \dot{\mathbf{a}} - \mathbf{\Omega} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{\Omega} = \mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F} \qquad (135) $$

Because the right hand member of Eq. (135) only depends of the rate of strain tensor that is objective, we can conclude on the objectivity of the equation governing the evolution of the second moment $\mathbf{a}$ of the orientation distribution $\psi(\mathbf{p})$.

Moreover, as $\mathbf{p} \cdot \mathbf{p} = 1$ and the orientation distribution function $\psi(\mathbf{p})$ satisfies the normality condition, it is immediate to prove that $\mathbf{a}$ has a unit trace, i.e. $\mathrm{Tr}(\mathbf{a}) = 1$. In order to ensure it, the terms involved in the equation of $\dot{\mathbf{a}}$ must be traceless. It is easy to verify that this is the case in Eq. (133).
### 5.3.2 On the absence of closure relations

Because the cluster is animated of a rigid motion, its kinematics is fully described by the second order moment **a** of the orientations distribution related to the rods that constitute it. At each time the orientation tensor **a** results a simple rotation of the initial one, and therefore, only the second order moment is involved, and no closure relations are needed for describing the microstructure evolution.

### 5.3.3 Numerical results

In what follows we consider a rod oriented initially along the direction $\mathbf{p}^0 = (\cos \theta^0, \sin \theta^0)^T$ immersed in a simple shear flow characterized by the homogeneous velocity field $\mathbf{v}^T = (\dot{\gamma} \cdot y, 0)$, with $\dot{\gamma} = -1$. At time $t$ the orientation of the rod will be given by $\theta$ (associated to the unit vector $\mathbf{p} = (\cos \theta, \sin \theta)^T$) and its rotational velocity by $\dot{\theta}$.

When the rod is assumed free of entanglements its kinematics is governed by the Jeffery's equation:

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (136) $$

Introducing the expressions of **p**

$$ \mathbf{p} = \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix} \quad (137) $$

and $\dot{\mathbf{p}}$

$$ \dot{\mathbf{p}} = \begin{pmatrix} -\dot{\theta} \cdot \sin \theta \\ \dot{\theta} \cdot \cos \theta \end{pmatrix} \quad (138) $$

into Eq. (136), then multiplying the first equation by $-\sin \theta$, the second one by $\cos \theta$ and then adding both, it results the Jeffery's rotary velocity $\dot{\theta}^J$

$$ \dot{\theta}^J = -\mathbf{S}_1^J \cdot \sin \theta + \mathbf{S}_2^J \cdot \cos \theta \quad (139) $$

where $\mathbf{S}_i^J$, $i = 1, 2$, are the components of the vector $\mathbf{S}^J$

$$ \mathbf{S}^J = \nabla \mathbf{v} \cdot \mathbf{p} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \cdot \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix} \quad (140) $$

Integrating Eq. (139) we can compute the orientation evolution of the free rod $\theta^J(t)$. Now, this orientation could be compared to the one resulting for the same rod but now considered making part of a rigid cluster.

We consider at the initial time different orientation tensors **a** having $\mathbf{p}^0$ as one of its principal directions. The orientation tensors having $\mathbf{p}^0$ as one of its principal directions are defined by

$$ \mathbf{a} = \alpha \cdot \begin{pmatrix} \cos^2(\theta^0) & \sin(\theta^0) \cdot \cos(\theta^0) \\ \sin(\theta^0) \cdot \cos(\theta^0) & \sin^2(\theta^0) \end{pmatrix} + (1-\alpha) \cdot \begin{pmatrix} \sin^2(\theta^0) & -\sin(\theta^0) \cdot \cos(\theta^0) \\ -\sin(\theta^0) \cdot \cos(\theta^0) & \cos^2(\theta^0) \end{pmatrix} \quad (141) $$

If $\alpha = 1$ all the rods composing the rigid cluster are aligned along the direction $\mathbf{p}^0$ and the resulting orientation evolution will coincide with the one predicted by the Jeffery's equation. When it is not the case, the orientation will be calculated from
$$ \dot{\mathbf{p}} = \mathbf{W} \cdot \mathbf{p} \qquad (142) $$

with

$$ \mathbf{W} = \mathbf{\Omega} + \mathbf{F} \qquad (143) $$

$$ \mathbf{F} = \begin{pmatrix} 0 & -f \\ f & 0 \end{pmatrix} \qquad (144) $$

and

$$ f = -\mathbf{D}_{12} - 2 \cdot \mathbf{D}_{11} \cdot \mathbf{a}_{12} + 2 \cdot \mathbf{D}_{12} \cdot \mathbf{a}_{11} \qquad (145) $$

Finally, the time evolution of tensor **a** reads:

$$ \dot{\mathbf{a}} = \mathbf{W} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{W}^T \qquad (146) $$

By introducing the expressions of **p** and $\dot{\mathbf{p}}$ given by Eqs. (137) and (138) respectively, into Eq. (142), multiplying the first resulting equation by $-\sin\theta$ and the second one by $\cos\theta$ and then adding both, it results

$$ \dot{\theta}^R = -\mathbf{S}_1^R \cdot \sin\theta + \mathbf{S}_2^R \cdot \cos\theta \qquad (147) $$

where $\mathbf{S}_i^R$, $i = 1, 2$, are the components of the vector $\mathbf{S}^R$

$$ \mathbf{S}^R = \mathbf{W} \cdot \mathbf{p} = \begin{pmatrix} 0 & \frac{\dot{\gamma}}{2} - f \\ -\frac{\dot{\gamma}}{2} + f & 0 \end{pmatrix} \cdot \begin{pmatrix} \cos\theta \\ \sin\theta \end{pmatrix} \qquad (148) $$

Integrating Eq. (147) we can compute the orientation evolution $\theta^R(t)$ of the rod belonging to the rigid cluster.

Both orientations $\theta^J(t)$ and $\theta^R(t)$ will be compared in the following cases:

- $\alpha = 1$. In that case both evolutions coincide.

- $\alpha = 0.98$. The orientation evolutions are depicted in Fig. 8. It can be noticed that the cluster induces a faster rotation and also that the rod can scape from the trap at $\theta = \pi$ because the effect of the other fibers not aligned along the direction $\theta = \pi$. On the contrary the free rod evolving according to the Jeffery's equation cannot overpass the orientation $\theta = \pi$.

- $\alpha = 0.75$. The orientation evolutions are depicted in Fig. 9. It can be noticed that the cluster being more isotropic, the rods oriented in other directions contribute to the cluster rotation that becomes the more and more uniform, and the traps located at $\theta = n \cdot \pi$, with $n = 1, 2, \dots$, are less significant.

- $\alpha = 0.5$. The orientation evolutions are depicted in Fig. 10. In this case because the isotropy of the cluster, the cluster rotary velocity is constant.
<figure><img src="image_9.png" /><figcaption>Fig. 8 Comparing $\theta^J(t)$ and $\theta^R(t)$ for a cluster characterized by $\alpha = 0.98$ in Eq. (141)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 Comparing $\theta^J(t)$ and $\theta^R(t)$ for a cluster characterized by $\alpha = 0.75$ in Eq. (141)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10 Comparing $\theta^J(t)$ and $\theta^R(t)$ for a cluster characterized by $\alpha = 0.5$ in Eq. (141)</figcaption></figure>
## 5.4 Kinetic theory description of a suspension of rigid clusters

The fact that rigid clusters rotate with a velocity that only depends on the second moment $\mathbf{a}$ of its orientation distribution function $\psi(\mathbf{p})$ suggests a microstructure description of suspensions of rigid clusters based on a distribution function $\Psi(\mathbf{x}, t, \mathbf{a})$ given the fraction of clusters that at position $\mathbf{x}$ and time $t$ have a configuration described by $\mathbf{a}$ related to the orientation of the rods composing it and that was noted by $\psi(\mathbf{p})$.

Function $\Psi(\mathbf{a})$ is subjected locally (at each position $\mathbf{x}$ and time $t$) to the normality condition:

$$ \int_A \Psi(\mathbf{a}) \, d\mathbf{a} = 1 \qquad (149) $$

The balance related to $\Psi$ writes:

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \Psi) + \frac{\partial}{\partial \mathbf{a}} (\dot{\mathbf{a}} \cdot \Psi) = 0 \qquad (150) $$

where

$$ \dot{\mathbf{a}} = \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F} \qquad (151) $$

with

$$ \mathbf{F} = \begin{pmatrix} 0 & -f \\ f & 0 \end{pmatrix} \qquad (152) $$

and $f$ given by

$$ f(\mathbf{D}, \mathbf{a}) = -\mathbf{D}_{12} - 2 \cdot \mathbf{D}_{11} \cdot \mathbf{a}_{12} + 2 \cdot \mathbf{D}_{12} \cdot \mathbf{a}_{11} \qquad (153) $$

### 5.4.1 Numerical results

In what follows we consider a suspension of rigid mono-disperse clusters in a simple shear flow characterized by the homogeneous velocity field $\mathbf{v}^T = (\dot{\gamma} \cdot \mathbf{y}, 0)$, with $\dot{\gamma} = 1$. The suspension microstructure is described by the distribution function $\Psi(\mathbf{x}, t, \mathbf{a})$. For the sake of simplicity we assume a homogeneous flow such that the microstructure does not depend of the physical space coordinates $\mathbf{x}$, being expressed by $\Psi(t, \mathbf{a})$. Moreover in 2D and because the symmetry of $\mathbf{a}$ and its unit trace, the microstructure will be described from only two components of the second order moment: in what follows we condier $\mathbf{a}_{11}$ and $\mathbf{a}_{12}$.

From

$$ \dot{\mathbf{a}} = \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F} \qquad (154) $$

we obtain after some trivial algebraic manipulations

$$ \dot{\mathbf{a}} = \begin{pmatrix} 2 \cdot \dot{\gamma} \cdot (1 - \mathbf{a}_{11}) \cdot \mathbf{a}_{12} & \dot{\gamma} \cdot (1 - 3 \cdot \mathbf{a}_{11} + 2 \cdot \mathbf{a}_{11}^2) \\ \dot{\gamma} \cdot (1 - 3 \cdot \mathbf{a}_{11} + 2 \cdot \mathbf{a}_{11}^2) & -2 \cdot \dot{\gamma} \cdot (1 - \mathbf{a}_{11}) \cdot \mathbf{a}_{12} \end{pmatrix} \qquad (155) $$

from which we can ensure that the evolution of the orientation tensor $\mathbf{a}$ will preserve the symmetry and its unit trace.

Now, the orientation distribution balance expressed by the Fokker-Planck equation (150) writes:

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{a}_{11}} (\dot{\mathbf{a}}_{11} \cdot \Psi) + \frac{\partial}{\partial \mathbf{a}_{12}} (\dot{\mathbf{a}}_{12} \cdot \Psi) = 0 \qquad (156) $$

that can be rewritten as:
$$ \frac{\partial \Psi}{\partial t} + \left( \frac{\partial \dot{\mathbf{a}}_{11}}{\partial \mathbf{a}_{11}} + \frac{\partial \dot{\mathbf{a}}_{12}}{\partial \mathbf{a}_{12}} \right) \cdot \Psi + \dot{\mathbf{a}}_{11} \cdot \frac{\partial \Psi}{\partial \mathbf{a}_{11}} + \dot{\mathbf{a}}_{12} \cdot \frac{\partial \Psi}{\partial \mathbf{a}_{12}} = 0 \quad (157) $$

where taking into account the expression of the components of $\dot{\mathbf{a}}$ given by (155) it results

$$ \frac{\partial \dot{\mathbf{a}}_{11}}{\partial \mathbf{a}_{11}} = -2 \cdot \dot{\gamma} \cdot \mathbf{a}_{12} \quad (158) $$

and

$$ \frac{\partial \dot{\mathbf{a}}_{12}}{\partial \mathbf{a}_{12}} = 0 \quad (159) $$

The resulting balance equation (157) must be solved in the conformation domain $\Omega = (0, 1) \times (-\frac{1}{2}, \frac{1}{2})$, i.e. $(\mathbf{a}_{11}, \mathbf{a}_{12}) \in \Omega$, but prior to apply any discretization technique we should define the boundary conditions.

For the discussion that follows we explicit the expression of the second order moment $\mathbf{a}$:

$$ \mathbf{a} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \psi(\mathbf{p}) \, d\mathbf{p} \quad (160) $$

and taking into account the expression of $\mathbf{p}$, $\mathbf{p}^T = (\cos\theta, \sin\theta)$, the components of $\mathbf{a}$ result:

$$ \mathbf{a} = \begin{pmatrix} \int_{0}^{2\pi} \cos^2 \theta \cdot \psi(\theta) \, d\theta & \int_{0}^{2\pi} \cos \theta \cdot \sin \theta \cdot \psi(\theta) \, d\theta \\ \int_{0}^{2\pi} \cos \theta \cdot \sin \theta \cdot \psi(\theta) \, d\theta & \int_{0}^{2\pi} \sin^2 \theta \cdot \psi(\theta) \, d\theta \end{pmatrix} \quad (161) $$

from which we can conclude that:

$$ \left\{ \begin{array}{l} 0 \le \mathbf{a}_{11} \le 1 \\ -\frac{1}{2} \le \mathbf{a}_{12} \le \frac{1}{2} \end{array} \right. \quad (162) $$

We analyze the limit cases:

- When $\mathbf{a}_{11} = 1$ the distribution function results necessarily $\psi(\theta) = \delta(\theta - 0)$, i.e. all the fibers align along the x-direction, and then $\mathbf{a}_{12} = 0$. Thus $\Psi(\mathbf{a}_{11} = 1, \mathbf{a}_{12} \neq 0) = 0$.

- When $\mathbf{a}_{11} = 0$ the distribution function results necessarily $\psi(\theta) = \delta(\theta - \frac{\pi}{2})$, i.e. all the fibers align along the y-direction, and then $\mathbf{a}_{12} = 0$. Thus $\Psi(\mathbf{a}_{11} = 0, \mathbf{a}_{12} \neq 0) = 0$.

- When $\mathbf{a}_{12} = \frac{1}{2}$ the distribution function results necessarily $\psi(\theta) = \delta(\theta - \frac{\pi}{4})$, i.e. all the fibers align along the $\frac{\pi}{4}$-direction, and then $\mathbf{a}_{11} = \frac{1}{2}$. Thus $\Psi(\mathbf{a}_{11} \neq \frac{1}{2}, \mathbf{a}_{12} = \frac{1}{2}) = 0$.

- When $\mathbf{a}_{12} = -\frac{1}{2}$ the distribution function results necessarily $\psi(\theta) = \delta(\theta + \frac{\pi}{4})$, i.e. all the fibers align along the $-\frac{\pi}{4}$-direction, and then $\mathbf{a}_{11} = \frac{1}{2}$. Thus $\Psi(\mathbf{a}_{11} \neq \frac{1}{2}, \mathbf{a}_{12} = -\frac{1}{2}) = 0$.

As there are 4 configurations on the domain boundary where the distribution function does not vanish necessarily, we should prove that in any case we can ensure the conservation balance. For this purpose it suffices to prove that the conformational velocity normal to the domain boundary at those points vanishes. Taking into account the expression of the conformational velocities given in Eq. (155):
<figure><img src="image_12.png" /><figcaption>Fig. 11 Velocity $\dot{\mathbf{a}}_{11}$ in $\Omega$</figcaption></figure>

$$ \left\{ \begin{aligned} \dot{\mathbf{a}}_{11} &= 2 \cdot \gamma \cdot (1 - \mathbf{a}_{11}) \cdot \mathbf{a}_{12} \\ \dot{\mathbf{a}}_{12} &= \gamma \cdot (1 - 3 \cdot \mathbf{a}_{11} + 2 \cdot \mathbf{a}_{11}^2) \end{aligned} \right. \qquad (163) $$

it results:

- If $\mathbf{a}_{11} = 0$ and $\mathbf{a}_{12} = 0$ then $\dot{\mathbf{a}}_{11} = 0$.

- If $\mathbf{a}_{11} = 1$ and $\mathbf{a}_{12} = 0$ then $\dot{\mathbf{a}}_{11} = 0$.

- If $\mathbf{a}_{11} = \frac{1}{2}$ and $\mathbf{a}_{12} = \frac{1}{2}$ then $\dot{\mathbf{a}}_{12} = 0$.

- If $\mathbf{a}_{11} = \frac{1}{2}$ and $\mathbf{a}_{12} = -\frac{1}{2}$ then $\dot{\mathbf{a}}_{12} = 0$.

that ensures the conservation balance.

The isotropic orientation distribution is defined by the center of the conformation domain $\Omega$, i.e. $\mathbf{a}_{11} = 0.5$ and $\mathbf{a}_{12} = 0$. The velocities $\dot{\mathbf{a}}_{11}$ and $\dot{\mathbf{a}}_{12}$ in the domain $\Omega$ given by Eq. (163) are depicted in Figs. 11 and 12 respectively, whereas the velocities directions are sketched in Fig. 13.

Because the balance equation involves pure advection appropriate stabilizations must be considered for integrating it. For example when using the finite differences, the first order derivatives of the distribution function must be up-winded, i.e. at node $(i, j)$ we should consider:

- If $(i, j)$ is located in $\mathbf{a}_{12} \ge 0$ where $\dot{\mathbf{a}}_{11} \ge 0$, then

$$ \left. \frac{\partial \Psi}{\partial \mathbf{a}_{11}} \right|_{i,j} \approx \frac{\Psi_{i,j} - \psi_{i-1,j}}{h_{11}} \qquad (164) $$

where $h_{11}$ is the distance between nodes in the $\mathbf{a}_{11}$-direction. If $(i, j)$ is located in $\mathbf{a}_{12} < 0$ where $\dot{\mathbf{a}}_{11} < 0$, then

$$ \left. \frac{\partial \Psi}{\partial \mathbf{a}_{11}} \right|_{i,j} \approx \frac{\Psi_{i+1,j} - \psi_{i,j}}{h_{11}} \qquad (165) $$
<figure><img src="image_13.png" /><figcaption>Fig. 12 Velocity $\dot{a}_{12}$ in $\Omega$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13 Directions of the conformational velocity</figcaption></figure>

- Analogously, when considering the velocity $\dot{a}_{12}$, if $(i, j)$ is located in $a_{11} \ge \frac{1}{2}$ where $\dot{a}_{12} \le 0$, then

$$ \left. \frac{\partial \Psi}{\partial \mathbf{a}_{12}} \right|_{i,j} \approx \frac{\Psi_{i,j+1} - \Psi_{i,j}}{h_{12}} \qquad (166) $$
<figure><img src="image_15.png" /><figcaption>Fig. 14 Time evolution of $\Psi(t, \mathbf{a}_{11}, \mathbf{a}_{12})$: (top-left) $\Psi(t = 0, \mathbf{a}_{11}, \mathbf{a}_{12})$; (top-right) $\Psi(t = 1, \mathbf{a}_{11}, \mathbf{a}_{12})$; (middle-left) $\Psi(t = 2, \mathbf{a}_{11}, \mathbf{a}_{12})$; (middle-right) $\Psi(t = 3, \mathbf{a}_{11}, \mathbf{a}_{12})$; (down-left) $\Psi(t = 4, \mathbf{a}_{11}, \mathbf{a}_{12})$; (down-right) $\Psi(t = 6, \mathbf{a}_{11}, \mathbf{a}_{12})$</figcaption></figure>

where $h_{12}$ is the distance between nodes in the $\mathbf{a}_{12}$-direction. If $(i, j)$ is located in $\mathbf{a}_{11} < \frac{1}{2}$ where $\dot{\mathbf{a}}_{12} > 0$, then

$$ \left. \frac{\partial \Psi}{\partial \mathbf{a}_{12}} \right|_{i,j} \approx \frac{\Psi_{i,j} - \psi_{i,j-1}}{h_{12}} \qquad (167) $$

We discretized the balance equation by using an implicit upwind finite differences scheme. However, despite its simplicity, the numerical scheme exhibited too high numerical diffusion for very localized initial orientation distributions and for long times. Fig. 14 depicts a sequence of snapshots of the time evolution of the distribution function $\Psi(t, \mathbf{a}_{11}, \mathbf{a}_{12})$.

To improve the accuracy avoiding artificial numerical diffusion, we decided to implement a discretization of the balance equation based on the method of particles [18] [19].

This technique consists in approximating the initial distribution $\Psi(t = 0, \mathbf{a})$ from $\mathcal{D}$ Dirac's masses:
$$ \Psi(t=0, \mathbf{a}) = \sum_{i=1}^{i=D} \alpha_i \cdot \delta(\mathbf{a} - \mathbf{a}_i^0) \qquad (168) $$

Now, the time evolution is computed by integrating

$$ \mathbf{a}_i(t) = \mathbf{a}_i^0 + \int_{\tau=0}^{\tau=t} \dot{\mathbf{a}}_i(\mathbf{a}_i(\tau)) \, d\tau \qquad (169) $$

and the orientation distribution at time t is reconstructed from

$$ \Psi(t, \mathbf{a}) = \sum_{i=1}^{i=D} \alpha_i \cdot \delta(\mathbf{a} - \mathbf{a}_i(t)) \qquad (170) $$

This technique, that in the incompressible case corresponds to an integration by characteristics, can be extended for considering diffusion terms as described in [5] by smoothing the Dirac's masses. However, as the equation that we are here addressing only involves pure advection, the standard method of particles runs perfectly. It can be noticed that this method is the equivalent to a direct integration of the Jeffery's equation when applied to the orientation distribution related to suspensions of ellipsoidal particles described by $\mathbf{p}$ or to an integration by characteristics when applied to the equation governing the evolution of the second order moment $\mathbf{a}$.

We first validate the integration scheme by considering an isotropic mono-disperse suspension of rigid clusters. In that case the initial distribution can be approximated by a single Dirac's mass concentrated at $\mathbf{a}^0 = (\mathbf{a}_{11}^0 = \frac{1}{2}, \mathbf{a}_{12}^0 = 0)$. In that case:

$$ \Psi(t=0, \mathbf{a}) = \delta(\mathbf{a} - \mathbf{a}^0) \qquad (171) $$

and because according to Eq. (163) the conformational velocity vanishes in the case of isotropic microstructures, then we obtain:

$$ \Psi(t, \mathbf{a}) = \delta(\mathbf{a} - \mathbf{a}^0), \quad \forall t \qquad (172) $$

When considering a mono-disperse microstructure characterized by $\mathbf{a}_{11}^0 > \frac{1}{2}$ and $\mathbf{a}_{12}^0 = 0$, the orientation tensor $\mathbf{a}$ evolves because the rotation of the aggregate. Some time later, a quarter of the cluster rotation period $T$, $t = \frac{T}{4}$, the cluster has rotated an angle $\frac{\pi}{2}$ and the orientation becomes:

$$ \mathbf{a} \left( t = \frac{T}{4} \right) = \begin{pmatrix} 1 - \mathbf{a}_{11}^0 & 0 \\ 0 & \mathbf{a}_{11}^0 \end{pmatrix} \qquad (173) $$

Obviously, after a half of a complete turn, $t = \frac{T}{2}$, the microstructure come back to the initial state, i.e. $\mathbf{a}(t = \frac{T}{2}) = \mathbf{a}^0$. Then $\mathbf{a}(t = n \cdot \frac{T}{2}) = \mathbf{a}^0$, with $n$ an arbitrary integer ($n \in \mathcal{N}$). Fig. 15 depicts the periodic trajectory in the space of configurations $(\mathbf{a}_{11}, \mathbf{a}_{12})$.

Finally we are considering an elongational flow characterized by the incompressible velocity field $\mathbf{v}^T = (\epsilon \cdot x, -\epsilon \cdot y)$.

From the expressions of $\Omega$ and $\mathbf{F}$ we can calculate $\mathbf{W}$ and from it the expression of $\dot{\mathbf{a}}$ that in the present case writes:

$$ \dot{\mathbf{a}} = \begin{pmatrix} 4 \cdot \epsilon \cdot \mathbf{a}_{12}^2 & 2 \cdot \epsilon \cdot \mathbf{a}_{12} \cdot (\mathbf{a}_{22} - \mathbf{a}_{11}) \\ 2 \cdot \epsilon \cdot \mathbf{a}_{12} \cdot (\mathbf{a}_{22} - \mathbf{a}_{11}) & -4 \cdot \epsilon \cdot \mathbf{a}_{12}^2 \end{pmatrix} \qquad (174) $$
<figure><img src="image_16.png" /><figcaption>Fig. 15 Trajectory in the conformational space related to $a_{11}^0 = 0.8$ and $a_{12}^0 = 0$.</figcaption></figure>

Thus, as the tensor describing the initial microstructure $a^0$ is symmetric and has a unit trace, its time evolution (174) will preserve both, the symmetry and the unit trace.

To verify the model we considered the flow characterized by $\epsilon = 1$ and an initial microstructure given by:

$$ \tilde{\mathbf{a}} = \begin{pmatrix} 0.8 & 0 \\ 0 & 0.2 \end{pmatrix} \qquad (175) $$

rotated of an angle $\theta = \pi/4$, i.e. $a^0 = \mathbf{R}^T \cdot \tilde{\mathbf{a}} \cdot \mathbf{R}$ with

$$ \mathbf{R} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{pmatrix} \qquad (176) $$

and $\theta = \pi/4$.

Thus, the initial microstructure reads:

$$ a^0 = \begin{pmatrix} 0.5 & 0.3 \\ 0.3 & 0.5 \end{pmatrix} \qquad (177) $$

It is expected that the elongation orients the rigid cluster in order to align its major axis along the elongation direction (x-direction) in order to obtain the steady solution

$$ a(t \to \infty) = \begin{pmatrix} 0.8 & 0 \\ 0 & 0.2 \end{pmatrix} \qquad (178) $$

where the cluster rotates without experiencing deformation because its rigid nature.

The trajectory $a(t)$ that the microstructure follows from $a^0$ until reaching the steady state $a(t \to \infty)$ is depicted in Fig. 16. This trajectory is defined by $a(t)$ computed by integrating
<figure><img src="image_17.png" /><figcaption>Fig. 16 Trajectory in the conformational space related to $a_{11}^0 = 0.5$ and $a_{12}^0 = 0.3$ when applying an elongational flow.</figcaption></figure>

$$ \mathbf{a}(t) = \mathbf{a}^0 + \int_{\tau=0}^{\tau=t} \dot{\mathbf{a}}(\mathbf{a}(\tau)) \, d\tau \qquad (179) $$

At each time, the microstructure is defined by the localized distribution function

$$ \Psi(t, \mathbf{a}) = \delta(\mathbf{a} - \mathbf{a}(t)) \qquad (180) $$

## 5.5 Virial stress

The contribution of a rigid cluster composed of N roods to the virial stress results:

$$ \tau = \xi \cdot L^2 \cdot \sum_{i=1}^{N} \mathbf{p}_i \otimes (\nabla \mathbf{v} \cdot \mathbf{p}_i - \dot{\mathbf{p}}_i) \qquad (181) $$

that taking into account the expression of $\dot{\mathbf{p}}_i$ reads:

$$ \tau = \xi \cdot L^2 \cdot \sum_{i=1}^{N} \mathbf{p}_i \otimes (\nabla \mathbf{v} \cdot \mathbf{p}_i - \mathbf{W} \cdot \mathbf{p}_i) \qquad (182) $$

or

$$ \tau = \xi \cdot L^2 \cdot \sum_{i=1}^{N} \left( (\mathbf{p}_i \otimes \mathbf{p}_i) \cdot (\nabla \mathbf{v})^T - (\mathbf{p}_i \otimes \mathbf{p}_i) \cdot \mathbf{W}^T \right) \qquad (183) $$

In the continuous mesoscopic framework Eq. (183) results:
$$ 
\begin{aligned}
\boldsymbol{\tau}^i &= 2 \cdot N_p \int_{S(0,1)} \left( (\mathbf{p} \otimes \mathbf{p}) \cdot (\nabla \mathbf{v})^T - (\mathbf{p} \otimes \mathbf{p}) \cdot \mathbf{W}^T \right) \cdot \psi(\mathbf{p}) \, d\mathbf{p} = \\
&= 2 \cdot N_p \cdot \left( \mathbf{a} \cdot (\nabla \mathbf{v})^T - \mathbf{a} \cdot \mathbf{W}^T \right) = 2 \cdot N_p \cdot \left( \mathbf{a} \cdot (\mathbf{D} + \mathbf{\Omega})^T - \mathbf{a} \cdot (\mathbf{\Omega} + \mathbf{F})^T \right) = \\
&= 2 \cdot N_p \cdot (\mathbf{a} \cdot \mathbf{D} + \mathbf{a} \cdot \mathbf{F})
\end{aligned}
\qquad (184) $$

The identity

$$ \mathbf{D} \cdot \mathbf{a} - \mathbf{F} \cdot \mathbf{a} = \mathbf{a} \cdot \mathbf{D} + \mathbf{a} \cdot \mathbf{F} \qquad (185) $$

ensures the stress tensor symmetry, i.e. $\boldsymbol{\tau} = \boldsymbol{\tau}^T$.

Now, by considering all the clusters that contribute to the local stress,

$$ \boldsymbol{\tau} = 2 \cdot \tilde{N}_p \cdot \int_{\mathcal{A}} (\mathbf{D} \cdot \mathbf{a} - \mathbf{F} \cdot \mathbf{a}) \cdot \Psi(\mathbf{a}) \, d\mathbf{a} \qquad (186) $$

In order to avoid spurious non-symmetries when solving the kinetic model by using stochastic simulations we could consider a symmetrized form of the virial stress:

$$ \boldsymbol{\tau} = \tilde{N}_p \cdot \int_{\mathcal{A}} (\mathbf{D} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{D} - \mathbf{F} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{F}) \cdot \Psi(\mathbf{a}) \, d\mathbf{a} \qquad (187) $$

where the fact that $\mathbf{F}^T = -\mathbf{F}$ has again being taken into account.

### 5.5.1 Model verification

#### Unidirectional suspension

There is a limit case in which the results already known should be found, the one consisting of rigid clusters, all them constituted of rods aligned in the same direction $\hat{\mathbf{p}}$. In that case, with $\psi(\mathbf{p}) = \delta(\mathbf{p} - \hat{\mathbf{p}})$, the classical Jeffery's theory predicts:

$$ \hat{\mathbf{a}} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \psi(\mathbf{p}) \, d\mathbf{p} = \hat{\mathbf{p}} \otimes \hat{\mathbf{p}} \qquad (188) $$

whose components write $\hat{\mathbf{a}}_{ij} = \hat{\mathbf{p}}_i \cdot \hat{\mathbf{p}}_j$.

Thus, we can explicitly calculate the expression $\mathbf{W} = \mathbf{\Omega} + \mathbf{F}$ and from it the one of $\hat{\mathbf{p}} = \mathbf{W} \cdot \hat{\mathbf{p}}$ and compare it with the one resulting from the Jeffery's expression $\hat{\mathbf{p}} = \nabla \mathbf{v} \cdot \hat{\mathbf{p}} - (\hat{\mathbf{p}}^T \cdot \mathbf{D} \cdot \hat{\mathbf{p}}) \cdot \hat{\mathbf{p}}$. As expected we obtain exactly the same expression.

#### Suspensions of discs

Other case of interest is the one consisting of rigid spheres, discs in 2D. As previously proved, the systems can be represented by a rigid cluster composed of two orthogonal rods having the same length.

In this case $\mathbf{F} = \mathbf{0}$ and then the virial stress reads:

$$ \boldsymbol{\tau} = 2 \cdot \gamma \cdot \mathbf{D} \qquad (189) $$

such that the suspension only experiences an increase of its viscosity $\tilde{\eta} = \eta + \gamma$ as predicted by the Einstein's formula.
# 6 Suspensions involving deformable clusters

In this section we are considering a more realistic scenario, that consisting in deformable clusters.

## 6.1 Microscopic kinematics

We assume two kind of applied forces, one due to the fluid-rod friction once more modeled from

$$ \mathbf{F}_i^H = \xi \cdot (\nabla \mathbf{v} \cdot \mathbf{p}_i \cdot L - \dot{\mathbf{p}}_i \cdot L) \qquad (190) $$

where the superscript "H" refers to its hydrodynamic nature.

On the other hand, we assume a force acting on each bead due to the rods entanglements. This force $\mathbf{F}_i^C$ is assumed scaling with the difference between the rigid motion velocity (the one that the bead would have if the cluster would be rigid) and the real one:

$$ \mathbf{F}_i^C = \mu \cdot (\mathbf{W} \cdot \mathbf{p}_i \cdot L - \dot{\mathbf{p}}_i \cdot L) \qquad (191) $$

By adding both forces it results:

$$ \mathbf{F}_i = \mathbf{F}_i^H + \mathbf{F}_i^C = L \cdot ((\mu \cdot \mathbf{W} + \xi \cdot \nabla \mathbf{v}) \cdot \mathbf{p}_i - (\xi + \mu) \cdot \dot{\mathbf{p}}_i) \qquad (192) $$

that can be rewritten as:

$$ \mathbf{F}_i = L \cdot (\xi + \mu) \cdot \left( \frac{\mu \cdot \mathbf{W} + \xi \cdot \nabla \mathbf{v}}{\xi + \mu} \cdot \mathbf{p}_i - \dot{\mathbf{p}}_i \right) \qquad (193) $$

Now, if we define the equivalent traceless gradient $\mathbf{G}$

$$ \mathbf{G} = \frac{\mu \cdot \mathbf{W} + \xi \cdot \nabla \mathbf{v}}{\xi + \mu} \qquad (194) $$

we can generalize all the results obtained for suspensions of dilute rods, by considering an equivalent friction coefficient $\xi + \mu$ and an equivalent gradient of velocities $\mathbf{G}$. In particular the rod rotary velocity results:

$$ \dot{\mathbf{p}}_i = \mathbf{G} \cdot \mathbf{p}_i - (\mathbf{p}_i^T \cdot \mathbf{G} \cdot \mathbf{p}_i) \cdot \mathbf{p}_i \qquad (195) $$

The term $\mathbf{p}_i^T \cdot \mathbf{G} \cdot \mathbf{p}_i$ deserves some additional comment. Being both $\mathbf{W}$ and $\mathbf{\Omega}$ skew-symmetric it results: $\mathbf{p}_i^T \cdot \mathbf{W} \cdot \mathbf{p}_i = 0$ and the same in the case of considering $\mathbf{\Omega}$, implying $\mathbf{p}_i^T \cdot \mathbf{G} \cdot \mathbf{p}_i = \frac{\xi}{\xi+\mu} \cdot \mathbf{p}_i^T \cdot \mathbf{D} \cdot \mathbf{p}_i$ and then:

$$ \dot{\mathbf{p}}_i = \mathbf{G} \cdot \mathbf{p}_i - \frac{\xi}{\xi + \mu} (\mathbf{p}_i^T \cdot \mathbf{D} \cdot \mathbf{p}_i) \cdot \mathbf{p}_i \qquad (196) $$

that can be rewritten as

$$ \begin{aligned} \dot{\mathbf{p}}_i &= \frac{\mu}{\xi + \mu} \cdot \mathbf{W} \cdot \mathbf{p}_i + \frac{\xi}{\xi + \mu} \cdot \nabla \mathbf{v} \cdot \mathbf{p}_i - \frac{\xi}{\xi + \mu} (\mathbf{p}_i^T \cdot \mathbf{D} \cdot \mathbf{p}_i) \cdot \mathbf{p}_i = \\ &= \frac{\xi}{\xi + \mu} \cdot \dot{\mathbf{p}}_i^J + \frac{\mu}{\xi + \mu} \cdot \mathbf{W} \cdot \mathbf{p}_i = \frac{\xi}{\xi + \mu} \cdot \dot{\mathbf{p}}_i^J + \frac{\mu}{\xi + \mu} \cdot \dot{\mathbf{p}}_i^R \end{aligned} \qquad (197) $$
where $\dot{\mathbf{p}}_i^J$ is the Jeffery's dilute contribution and $\dot{\mathbf{p}}_i^R$ the one coming from the rigid cluster contribution.

We can notice that when $\xi \gg \mu$ the cluster does not resist the flow induced rod orientation, and in that case $\dot{\mathbf{p}}_i \approx \dot{\mathbf{p}}_i^J$. In the opposite case, $\mu \gg \xi$, the cluster is too rigid and the bead adopts the velocity dictated by the rigid cluster kinematics $\dot{\mathbf{p}}_i \approx \dot{\mathbf{p}}_i^R$.

## 6.2 Mesoscopic description

If we assume the orientation of the rods of a deformable cluster given by the orientation distribution $\psi(\mathbf{p})$, coarser descriptions will need the use of some of its moments as well as the equations governing their time evolution.

We start considering again the usual second order moment:

$$ \mathbf{a} = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \psi(\mathbf{p}) \, d\mathbf{p} \qquad (198) $$

whose time derivative reads

$$ \dot{\mathbf{a}} = \int_{S(0,1)} (\dot{\mathbf{p}} \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}) \cdot \psi(\mathbf{p}) \, d\mathbf{p} \qquad (199) $$

where as just deduced the microscopic velocity is given by:

$$ \dot{\mathbf{p}} = \frac{\xi}{\xi + \mu} \cdot \dot{\mathbf{p}}^J + \frac{\mu}{\xi + \mu} \cdot \dot{\mathbf{p}}^R \qquad (200) $$

that introduced into the expression of the second order moment time derivative results

$$ \dot{\mathbf{a}} = \frac{\xi}{\xi + \mu} \cdot \dot{\mathbf{a}}^J + \frac{\mu}{\xi + \mu} \cdot \dot{\mathbf{a}}^R \qquad (201) $$

with $\dot{\mathbf{a}}^J$ and $\dot{\mathbf{a}}^R$ given by Eqs. (32) and (111) respectively:

$$ \left\{ \begin{aligned} \dot{\mathbf{a}}^J &= \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \mathbf{D} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{D} - 2 \cdot \mathbf{A} : \mathbf{D} \\ \dot{\mathbf{a}}^R &= \mathbf{W} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{W}^T = \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F} \end{aligned} \right. \qquad (202) $$

Introducing Eq. (202) into Eq. (201) and grouping the different terms it results:

$$ \dot{\mathbf{a}} = \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \frac{\xi}{\xi + \mu} \cdot (\mathbf{D} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{D} - 2 \cdot \mathbf{A} : \mathbf{D}) + \frac{\mu}{\xi + \mu} \cdot (\mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F}) \qquad (203) $$

that can be rewritten as:

$$ \begin{aligned} \dot{\mathbf{a}} - \mathbf{\Omega} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{\Omega} - \frac{\xi}{\xi + \mu} \cdot (\mathbf{D} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{D}) &= \\ &= -\frac{\xi}{\xi + \mu} \cdot 2 \cdot \mathbf{A} : \mathbf{D} + \frac{\mu}{\xi + \mu} \cdot (\mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F}) \end{aligned} \qquad (204) $$

where the left hand member defines the so called Johnson & Segalman objective derivative:

$$ \frac{\delta_c(\bullet)}{\delta t} = (\dot{\bullet}) - \mathbf{\Omega} \cdot (\bullet) + (\bullet) \cdot \mathbf{\Omega} - c \cdot (\mathbf{D} \cdot (\bullet) + (\bullet) \cdot \mathbf{D}) \qquad (205) $$

for $c = \frac{\xi}{\xi+\mu}$.

The two limit cases are:
- When flow induced forces are preponderant, $\xi \gg \mu$ and $c \approx 1$. In that case de Johnson & Segalman derivative approaches de Oldroyd's one that appears in the standard Jeffery's based rods suspension modeling.

- When cluster effects are preponderant, $\mu \gg \xi$ and then the Johnson & Segalman derivative approaches de Jauman's one characterized by $c = 0$.

In any case, we just proved that the model just derived that incorporates drag and network effects is fully objective.

### 6.2.1 Numerical results

In what follows we consider again a rod oriented initially along the direction $\mathbf{p}^0 = (\cos \theta^0, \sin \theta^0)^T$ immersed in a simple shear flow characterized by the homogeneous velocity field $\mathbf{v}^T = (\dot{\gamma} \cdot \mathbf{y}, 0)$, with $\dot{\gamma} = -1$. At time $t$ the orientation of the rod will be given by $\theta$ (associated to the unit vector $\mathbf{p} = (\cos \theta, \sin \theta)^T$) and its rotational velocity by $\dot{\theta}$.

Again, when the rod is assumed free of entanglements its kinematics is governed by the Jeffery's equation:

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (206) $$

that as previously deduced results in:

$$ \dot{\theta}^J = -\mathbf{S}_1^J \cdot \sin \theta + \mathbf{S}_2^J \cdot \cos \theta \quad (207) $$

whose integration leads the orientation evolution of the free rod $\theta^J(t)$. Now, this orientation will be compared to the one resulting for the same rod but now considered making part of a deformable cluster.

We consider again at the initial time different orientation tensors $\mathbf{a}$ having $\mathbf{p}^0$ as one of its principal directions, given by

$$ \mathbf{a} = \alpha \cdot \begin{pmatrix} \cos^2(\theta^0) & \sin(\theta^0) \cdot \cos(\theta^0) \\ \sin(\theta^0) \cdot \cos(\theta^0) & \sin^2(\theta^0) \end{pmatrix} + (1-\alpha) \begin{pmatrix} \sin^2(\theta^0) & -\sin(\theta^0) \cdot \cos(\theta^0) \\ -\sin(\theta^0) \cdot \cos(\theta^0) & \cos^2(\theta^0) \end{pmatrix} \quad (208) $$

If $\alpha = 1$ all the rods composing the cluster are aligned along the direction $\mathbf{p}^0$ and again the resulting orientation evolution will coincide with the one predicted by the Jeffery's equation. When it is not the case, the orientation will be calculated from

$$ \dot{\mathbf{p}} = \mathbf{G} \cdot \mathbf{p} - \frac{\xi}{\xi + \mu} (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (209) $$

with

$$ \mathbf{G} = \frac{\xi \cdot \nabla \mathbf{v} + \mu \cdot \mathbf{W}}{\xi + \mu} \quad (210) $$

$$ \mathbf{W} = \mathbf{\Omega} + \mathbf{F} \quad (211) $$

$$ \mathbf{F} = \begin{pmatrix} 0 & -f \\ f & 0 \end{pmatrix} \quad (212) $$
<figure><img src="image_18.png" /><figcaption>Fig. 17 Comparing $\theta^J(t)$ and $\theta^D(t)$ for a cluster characterized by $\xi = 0.5$, $\mu = 0.5$ and $\alpha = 0.75$</figcaption></figure>

and

$$f = -\mathbf{D}_{12} - 2 \cdot \mathbf{D}_{11} \cdot \mathbf{a}_{12} + 2 \cdot \mathbf{D}_{12} \cdot \mathbf{a}_{11} \quad (213)$$

Finally, the time evolution of tensor **a** reads:

$$\begin{aligned} \dot{\mathbf{a}} = \frac{\mu}{\xi + \mu} \cdot (\mathbf{W} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{W}^T) + \\ + \frac{\xi}{\xi + \mu} \cdot (\nabla \mathbf{v} \cdot \mathbf{a} + \mathbf{a} \cdot (\nabla \mathbf{v})^T - 2 \cdot (\mathbf{a} : \mathbf{D}) \cdot \mathbf{a}) \end{aligned} \quad (214)$$

where a quadratic closure relation has been considered.

By introducing the expressions of **p** and $\dot{\mathbf{p}}$ given by Eqs. (137) and (138) respectively, into Eq. (209), multiplying the first resulting equation by $-\sin\theta$ and the second one by $\cos\theta$ and then adding both, it results

$$\dot{\theta}^D = -\mathbf{S}_1^D \cdot \sin\theta + \mathbf{S}_2^D \cdot \cos\theta \quad (215)$$

where $\mathbf{S}_i^D$, $i = 1, 2$, are the components of the vector $\mathbf{S}^D$

$$\mathbf{S}^D = \mathbf{G} \cdot \mathbf{p} \quad (216)$$

Integrating Eq. (215) we can compute the orientation evolution of the rod belonging to the deformable cluster $\theta^D(t)$.

Both orientations $\theta^J(t)$ and $\theta^D(t)$ will be compared in the following cases:

- $\alpha = 1$. In that case both evolutions coincide for any value of $\xi$ and $\mu$.

- Figs. 17, 18 and 19 compare both orientation evolutions when $\xi = 0.5$, $\mu = 0.5$ and the rods constituting the initial cluster at time $t = 0$ describe an orientation tensor given by $\alpha = 0.75$, $\alpha = 0.5$ and $\alpha = 0.25$ respectively in Eq. (208).

- Figs. 20, 21 and 22 compare both orientation evolutions when $\xi = 0.25$, $\mu = 0.75$ and the rods constituting the initial cluster at time $t = 0$ describe an orientation tensor given by $\alpha = 0.75$, $\alpha = 0.5$ and $\alpha = 0.25$ respectively in Eq. (208).

From these results we can conclude that as expected the cluster effect increases the averaged rotary velocity.
<figure><img src="image_19.png" /><figcaption>**Fig. 18** Comparing $\theta^J(t)$ and $\theta^D(t)$ for a cluster characterized by $\xi = 0.5$, $\mu = 0.5$ and $\alpha = 0.5$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19** Comparing $\theta^J(t)$ and $\theta^D(t)$ for a cluster characterized by $\xi = 0.5$, $\mu = 0.5$ and $\alpha = 0.25$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 20** Comparing $\theta^J(t)$ and $\theta^D(t)$ for a cluster characterized by $\xi = 0.25$, $\mu = 0.75$ and $\alpha = 0.75$</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 21 Comparing $\theta^J(t)$ and $\theta^D(t)$ for a cluster characterized by $\xi = 0.25$, $\mu = 0.75$ and $\alpha = 0.5$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22 Comparing $\theta^J(t)$ and $\theta^D(t)$ for a cluster characterized by $\xi = 0.25$, $\mu = 0.75$ and $\alpha = 0.25$</figcaption></figure>

## 6.3 Kinetic description

In the present case, by assuming an appropriate closure relation for describing the fourth moment $\mathbf{A}$ of the orientation distribution, we can consider again a kinetic theory description based on the distribution function $\Psi(\mathbf{x}, t, \mathbf{a})$ given the fraction of clusters that at position $\mathbf{x}$ and time $t$ have a configuration described by $\mathbf{a}$ related to the orientation of the rods composing it and that was noted by $\psi(\mathbf{p})$.

Function $\Psi(\mathbf{a})$ is subjected locally (at each position $\mathbf{x}$ and time $t$) to the normality condition:

$$ \int_{\mathcal{A}} \Psi(\mathbf{a}) \, d\mathbf{a} = 1 \qquad (217) $$

The balance related to $\Psi$ writes:

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \Psi) + \frac{\partial}{\partial \mathbf{a}} (\dot{\mathbf{a}} \cdot \Psi) = 0 \qquad (218) $$
<figure><img src="image_24.png" /><figcaption>Fig. 23 Trajectory in the conformational space related to $a_{11}^0 = 0.8$ and $a_{12} = 0$ for a deformable cluster.</figcaption></figure>

where

$$ \dot{\mathbf{a}} = \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \frac{\xi}{\xi + \mu} \cdot (\mathbf{D} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{D} - 2 \cdot \mathbf{A} : \mathbf{D}) + \frac{\mu}{\xi + \mu} \cdot (\mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F}) \quad (219) $$

with

$$ \mathbf{F} = \begin{pmatrix} 0 & -f \\ f & 0 \end{pmatrix} \qquad (220) $$

and $f$ given by

$$ f(\mathbf{D}, \mathbf{a}) = -\mathbf{D}_{12} - 2 \cdot \mathbf{D}_{11} \cdot \mathbf{a}_{12} + 2 \cdot \mathbf{D}_{12} \cdot \mathbf{a}_{11} \quad (221) $$

### 6.3.1 Numerical results

In this section we are addressing the same numerical experiments that were considered in section 5.4.1 but now considering deformable clusters characterized by $\xi = 0.3$ and $\mu = 0.7$.

When considering a mono-disperse microstructure characterized by $a_{11}^0 > \frac{1}{2}$ and $a_{12}^0 = 0$ in a simple shear flow with $\dot{\gamma} = 1$, the orientation tensor $\mathbf{a}$ evolves because the rotation of the aggregate and its deformation. In these flow conditions given by the values of $\dot{\gamma}$, $\xi$ and $\mu$ the time evolution of the orientation tensor is again periodic in the conformation space. Because the flow induced cluster deformation the microstructure never reach the initial microstructure rotated of an angle $\frac{\pi}{2}$ as was the case when dealing with rigid clusters. Fig. 23 depicts the trajectory in the space of configurations $(\mathbf{a}_{11}, \mathbf{a}_{12})$.
Finally we are considering an elongational flow characterized by the incompressible velocity field $\mathbf{v}^T = (\epsilon \cdot x, -\epsilon \cdot y)$, with $\epsilon = 1$ and an initial microstructure given by

$$ \tilde{\mathbf{a}} = \begin{pmatrix} 0.8 & 0 \\ 0 & 0.2 \end{pmatrix} \qquad (222) $$

rotated of an angle $\theta = \frac{\pi}{4}$, i.e. $\mathbf{a}^0 = \mathbf{R}^T \cdot \tilde{\mathbf{a}} \cdot \mathbf{R}$ with

$$ \mathbf{R} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{pmatrix} \qquad (223) $$

and $\theta = \frac{\pi}{4}$. Thus, the initial microstructure reads:

$$ \mathbf{a}^0 = \begin{pmatrix} 0.5 & 0.3 \\ 0.3 & 0.5 \end{pmatrix} \qquad (224) $$

It is expected that the elongation orients the rigid cluster in order to align its major axis along the elongation direction (x-direction) but at the same time the elongation tends to extend the cluster to reach a full alignment of all the rods composing it in the direction of the applied extension, that results in the steady state orientation tensor:

$$ \mathbf{a}(t \to \infty) = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \qquad (225) $$

The trajectory $\mathbf{a}(t)$ that the microstructure follows from $\mathbf{a}^0$ until reaching the steady state $\mathbf{a}(t \to \infty)$ is depicted in Fig. 24. This trajectory is defined by $\mathbf{a}(t)$ computed by integrating

$$ \mathbf{a}(t) = \mathbf{a}^0 + \int_{\tau=0}^{\tau=t} \dot{\mathbf{a}}(\mathbf{a}(\tau)) \, d\tau \qquad (226) $$

At each time, the microstructure is defined by the localized distribution function

$$ \Psi(t, \mathbf{a}) = \delta(\mathbf{a} - \mathbf{a}(t)) \qquad (227) $$

Finally, in order to conclude on the qualitative and quantitative significance of the cluster effect, we compute the time evolution of a mono-disperse microstructure characterized by $\mathbf{a}_{11}^0 > \frac{1}{2}$ and $\mathbf{a}_{12}^0 = 0$ in a simple shear flow with $\dot{\gamma} = 1$ but neglecting the effects of the whole cluster on the orientation of each rod composing it, that is, by assuming $\mu = 0$ and $\xi = 1$. When starting from an initial microstructure characterized by $(\mathbf{a}_{11}^0 = 0.8, \mathbf{a}_{12}^0 = 0)$, the computed trajectory in the conformation space $(\mathbf{a}_{11}, \mathbf{a}_{12})$ is depicted in Fig. 25.

We can notice that the orientation of the rods composing the soft cluster evolves because the flow induced orientation, without any affect of the cluster because we set $\mu = 0$. When the orientation of rods reach the shear direction $\mathbf{p}^T = (1, 0)$ the rotation velocity vanishes and the rods remain at rest. Thus, the steady state microstructure is characterized by all the rods composing the cluster aligned in the shear direction, i.e. $\mathbf{a}_{11}(t \to \infty) = 1$ and $\mathbf{a}_{12}(t \to \infty) = 0$.
<figure><img src="image_25.png" /><figcaption>Fig. 24 Trajectory in the conformational space related to $a_{11}^0 = 0.5$ and $a_{12} = 0.3$ when applying an elongational flow.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 25 Trajectory in the conformational space related to $a_{11}^0 = 0.8$ and $a_{12} = 0$ in a simple shear flow when cluster effects are neglected.</figcaption></figure>
## 6.4 Virial stress

The virial stress associated with a deformable cluster is computed by considering the hydrodynamic forces acting on each bead $L \cdot \mathbf{p}$

$$ \mathbf{F}^H(\mathbf{p}) = \xi \cdot L \cdot (\nabla \mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}) \qquad (228) $$

from

$$ \tau_i = L \cdot \int_{S(0,1)} \mathbf{p} \otimes \mathbf{F}^H(\mathbf{p}) \cdot \psi_i(\mathbf{p}) \, d\mathbf{p} \qquad (229) $$

where the subindex "i" refer to the contribution of cluster i located around point x and time t where the stress is being computed.

If we introduce now into the expression of the hydrodynamic force the expression of the rod velocity given by Eq. (196)

$$ \dot{\mathbf{p}} = \mathbf{G} \cdot \mathbf{p} - \frac{\xi}{\xi + \mu} (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} \qquad (230) $$

it results

$$ \mathbf{F}^H(\mathbf{p}) = \xi \cdot L \cdot (\nabla \mathbf{v} \cdot \mathbf{p} - \mathbf{G} \cdot \mathbf{p} + \frac{\xi}{\xi + \mu} (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p}) \qquad (231) $$

Now, by making use of the relations:

$$ \left\{ \begin{array}{l} \mathbf{G} = \frac{\mu \cdot \mathbf{W} + \xi \cdot \nabla \mathbf{v}}{\xi + \mu} \\ \mathbf{W} = \mathbf{\Omega} + \mathbf{F} \\ \nabla \mathbf{v} = \mathbf{D} + \mathbf{\Omega} \end{array} \right. \qquad (232) $$

the expression of the hydrodynamic forces reduces to:

$$ \mathbf{F}^H(\mathbf{p}) = \xi \cdot L \cdot \left( \frac{\mu}{\xi + \mu} \cdot (\mathbf{D} \cdot \mathbf{p} - \mathbf{F} \cdot \mathbf{p}) + \frac{\xi}{\xi + \mu} (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p}) \cdot \mathbf{p} \right) \qquad (233) $$

from with the cluster virial becomes:

$$ \begin{aligned} \tau_i &= L \cdot \int_{S(0,1)} \mathbf{p} \otimes \mathbf{F}^H(\mathbf{p}) \cdot \psi_i(\mathbf{p}) \, d\mathbf{p} = \\ &= 2 \cdot N_p \cdot \left( \frac{\mu}{\xi + \mu} \cdot \boldsymbol{\tau}_i^R + \frac{\xi}{\xi + \mu} \cdot \boldsymbol{\tau}_i^J \right) \end{aligned} \qquad (234) $$

where the cluster and hydrodynamic contributions, $\boldsymbol{\tau}_i^R$ and $\boldsymbol{\tau}_i^J$, read:

$$ \boldsymbol{\tau}_i^R = \mathbf{D} \cdot \mathbf{a} - \mathbf{F} \cdot \mathbf{a} \qquad (235) $$

and

$$ \boldsymbol{\tau}_i^J = \mathbf{A} : \mathbf{D} \qquad (236) $$

respectively.

As previously proved, both terms are symmetric, ensuring the whole symmetry of the virial stress.

The two limits cases previously addressed, the Jeffery's regime is reached by assuming the cluster contribution negligible. In the opposite way, when cluster effects is dominant, the virial related to a rigid cluster is obtained.
## 7 Fist steps towards more concentrated regimes involving cluster-cluster interaction

When the concentration of aggregates is high enough, cluster interaction cannot be more neglected.

From one side, big aggregates break because the flow induced hydrodynamic forces. The destruction rate could be considered scaling with the equivalent rate of strain. We could assume in first approximation that a cluster breaks originating two smaller descendant clusters, having the same microstructure that the father cluster.

On the other side the coalescence is assumed also being produced from the interaction (collision) of two clusters for giving a bigger one. In general the coalescence rate is expected scaling with the equivalent rate of strain, because it depends on the cluster interactions that scales with the gradient of velocity. The resulting microstructure deserves some additional comment.

If cluster A composed of $N_A$ rods and having a microstructure described by the orientation distribution function $\psi_A(\mathbf{p})$ meets cluster B composed of $N_B$ rods and having a microstructure described by the orientation distribution function $\psi_B(\mathbf{p})$, the resulting normalized orientation distribution for the resulting cluster C reads:

$$ \psi_C(\mathbf{p}) = \frac{N_A \cdot \psi_A(\mathbf{p}) + N_B \cdot \psi_B(\mathbf{p})}{N_A + N_B} \qquad (237) $$

If clusters A and B were characterized by the second moment of their respective orientation distributions, $\mathbf{a}_A$ and $\mathbf{a}_B$ respectively, the moment for the resulting cluster C writes:

$$ \begin{aligned} \mathbf{a}_C(\mathbf{p}) &= \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \psi_C(\mathbf{p}) \, d\mathbf{p} = \\ &= \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \cdot \frac{N_A \cdot \psi_A(\mathbf{p}) + N_B \cdot \psi_B(\mathbf{p})}{N_A + N_B} \, d\mathbf{p} = \frac{N_A \cdot \mathbf{a}_A + N_B \cdot \mathbf{a}_B}{N_A + N_B} \end{aligned} \qquad (238) $$

that constitutes a simple mixing rule.

The kinetic modelling should integrate an additional coordinate quantifying the cluster size $S$ (related to the number of rods that it contains), allowing the mesoscopic description of polidisperse suspensions involving rigid or deformable clusters. Such microstructural description uses the distribution function $\Psi(\mathbf{x}, t, \mathbf{a}, S)$. When clusters involve too many rods, a continuous description of its size is more suitable than a discrete one based on the number of rods, because it allows defining coarse and fine microstructure descriptions depending on the level of discretization of the continuous coordinate $S$.

If we denote by $v_c$ and $v_d$ the cluster construction and destruction rates respectively, the associated terms in the distribution balance at position $\mathbf{x}$ and time $t$ will be denoted by $\dot{\Psi}_C(\mathbf{x}, t, \mathbf{a}, S)$ and $\dot{\Psi}_D(\mathbf{x}, t, \mathbf{a}, S)$ respectively. The resulting Fokker-Planck equation reads:

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \Psi) + \frac{\partial}{\partial \mathbf{a}} (\dot{\mathbf{a}} \cdot \Psi) = \dot{\Psi}_C(\mathbf{x}, t, \mathbf{a}, S) + \dot{\Psi}_D(\mathbf{x}, t, \mathbf{a}, S) \qquad (239) $$

where in the general case

$$ \dot{\mathbf{a}} = \mathbf{\Omega} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{\Omega} + \frac{\xi}{\xi + \mu} \cdot (\mathbf{D} \cdot \mathbf{a} + \mathbf{a} \cdot \mathbf{D} - 2 \cdot \mathbf{A} : \mathbf{D}) + \frac{\mu}{\xi + \mu} \cdot (\mathbf{F} \cdot \mathbf{a} - \mathbf{a} \cdot \mathbf{F}) \qquad (240) $$
with

$$ \mathbf{F} = \begin{pmatrix} 0 & -f \\ f & 0 \end{pmatrix} \qquad (241) $$

and $f$ given by

$$ f(\mathbf{D}, \mathbf{a}) = -\mathbf{D}_{12} - 2 \cdot \mathbf{D}_{11} \cdot \mathbf{a}_{12} + 2 \cdot \mathbf{D}_{12} \cdot \mathbf{a}_{11} \qquad (242) $$

The source terms related to the construction mechanism $\dot{\Psi}_C(\mathbf{x}, t, \mathbf{a})$ reads:

$$ \begin{aligned} \dot{\Psi}_C(\mathbf{x}, t, \mathbf{a}, S) = & \int_A \int_{S_{min}}^{S_{max}} v_C \cdot \Psi(\mathbf{x}, t, \mathbf{a}', r) \cdot \Psi(\mathbf{x}, t, \mathbf{a}'', S-r) \, dr \cdot d\mathbf{a} - \\ & - \int_A \int_{S_{min}}^{S_{max}-S} v_C \cdot \Psi(\mathbf{x}, t, \mathbf{a}, S) \cdot \Psi(\mathbf{x}, t, \mathbf{a}', r) \, dr \cdot d\mathbf{a} \end{aligned} \qquad (243) $$

where in order to ensure that the resulting microstructure generated from clusters having microstructures $\mathbf{a}'$ and $\mathbf{a}''$ is precisely $\mathbf{a}$, the microstructure $\mathbf{a}''$ must verify

$$ \frac{r \cdot \mathbf{a}' + (S-r) \cdot \mathbf{a}''}{S} = \mathbf{a} \qquad (244) $$

The first term in the right hand side of Eq. (243) express that two clusters of sizes $r$ and $S-r$ having appropriate microstructures meets with a probability depending on their fractions (expressed by the value of the corresponding distribution functions) to create a bigger cluster of size $S$ with a microstructure $\mathbf{a}$. However, clusters of size $S$ and microstructure $\mathbf{a}$ that meet any other cluster with compatible size $S_{min} \le r \le S_{max}-S$ and an arbitrary microstructure $\mathbf{a}'$ produce the decrease of population $(\mathbf{a}, S)$.

On the other hand, the source terms related to the destruction mechanism $\dot{\Psi}_D(\mathbf{x}, t, \mathbf{a})$ is given by:

$$ \begin{aligned} \dot{\Psi}_D(\mathbf{x}, t, \mathbf{a}, S) = & \int_S^{S_{max}} v_D \cdot \frac{1}{r} \cdot \Psi(\mathbf{x}, t, \mathbf{a}, r) \, dr - \\ & - v_D \cdot \Psi(\mathbf{x}, t, \mathbf{a}, S) \end{aligned} \qquad (245) $$

In this expression the first term accounts for the generation of clusters of size $S$ and microstructure $\mathbf{a}$ due to the destruction of bigger clusters with the same microstructure (we assumed in first approximation that the descendants have the same microstructure that the progenitor cluster). By assuming that all the sizes of the descendants are equiprobable, the probability that a cluster of size $r$ generates a cluster of size $z$ and other with the complementary size $r-z$ will be $\frac{1}{r}$ i.e. $\mathcal{P}(r; z, r-z) = \frac{1}{r}$, justifying the factor that appears in the first term of the right hand member of Eq. (245). The last term is trivial, because it is related to the simple destruction of cluster of size $S$ and microstructure $\mathbf{a}$.

It is easy to prove that both, the construction and destruction terms just proposed allow to ensure the conservation balance.

The main issue related to enriched kinetic theory descriptions lies in its curse of dimensionality, because generalized Fokker-Planck equations are defined in highly multidimensional spaces involving the physical coordinates (space and time) and a number of conformational coordinates describing the microstructure at the desired level of detail. Classically the curse of dimensionality was circumvented by different ways, being the most usual:
- the solution of the stochastic counterpart, however, it is well known that the control of the statistical noise and the variance reduction are the main tricky issues, making possible the calculation of the distribution moments but making difficult the reconstruction of the distribution itself;

- the derivation of macroscopic approaches, some times purely phenomenological, other times inspired or derived from finer descriptions, but even in the last case the use of closure approximations is generally compulsory;

- introducing a first discretization by substituting conformational coordinates by a series of populations. Thus, for example, in the polydisperse case we could introduce different populations of size $S_i$ and describe the microstructure of each one from its distribution function $\Psi_i(\mathbf{x}, t, \mathbf{a})$ whose evolution is given by the coupled Fokker-Planck equations:

$$ \frac{\partial \Psi_i}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\mathbf{v} \cdot \Psi_i) + \frac{\partial}{\partial \mathbf{a}} (\dot{\mathbf{a}} \cdot \Psi_i) = \sum_{j \neq i} (\Phi_{j \to i} - \Phi_{i \to j}) \quad (246) $$

where the source terms couple the different populations through the fluxes coming in or coming-out of each population $i$. Obviously, one could imagine to remove the second order moment as coordinate by considering additional populations characterized by their size $S_i$ and their microstructure $\mathbf{a}_k$ and described at the mesoscopic scale from the distribution functions $\Psi_{ik}(\mathbf{x}, t)$. This approach was considered in the kinetic description of systems of active particles, in order to circumvent the curse of dimensionality [14], however, as soon as the dimensionality increases the number of populations explodes;

- the use of separated representations in order to ensure that the complexity scales linearly with the model dimensionality. These approaches were widely employed in our former works, and because in our opinion it is the most plausible way to address such multidimensional models we revisits its foundations in the next section.

## 8 Efficient solvers of multidimensional partial differential equations

When we are addressing models in quantum chemistry (the finest description of atomic and molecular bonds responsible of the structure and mechanics of materials), the so called wavefunction $\Psi$ representing the distribution of an elementary particle is defined in the whole physical space, i.e. $\Psi(\mathbf{x}, t)$. Its evolution is governed by the Schrodinger equation (or its fully relativistic Dirac counterpart) that defines a standard 3D transient model. However, when the physical system contains $d$ particles, the evolution of the associated wavefunction $\Psi(\mathbf{x}_1, \dots, \mathbf{x}_d, t)$ is governed by a transient problem defined in a space involving $3 + 1 + 3 \cdot d$ dimensions [17] [43]. If one proceeds to the solution of such model by using a standard mesh based discretization technique, where $M$ nodes are used for discretizing each space $\mathbf{x}_i$, $\forall i \in [1, \dots, d]$, the resulting number of nodes reaches the astronomical value of $M^d$. With $M \approx 1000$ (a very coarse description in practice) and $d \approx 30$ (a very small atomic system) the numerical complexity results $10^{90}$. It is important to recall that this value corresponds to the presumed number of elementary particles in the universe. For this reason, the Nobel Prize R.B. Laughlin affirmed [42]:

> "... the theory of Everything is not even remotely a theory of everything. We know this equation (Laughlin refers to the Schrödinger equation) is correct because it has been
solved accurately for small number of particles -isolated atoms and small molecules- and found to agree in minute detail with experiments. However it cannot be solved accurately when the number of particles exceeds about 10. No computer existing, or that will ever exist, can break this barrier because it is a catastrophe of dimension. If the amount of computer memory required to represent the quantum wavefunction of one particle is $M$ then the amount required to represent the wavefunction of $d$ particles is $M^d$ ...

The curse of dimensionality is present in many other scales in the description of the structure and mechanics of materials. For example, when we consider a linear macromolecule, its conformation can be expressed by giving the position $\mathbf{r}_i$ of some points (beads) along its contour, or equivalently, by giving the vectors $\mathbf{q}_i$ joining two consecutive beads. Within the kinetic theory framework, molecular systems are described (see for example [15]) by a distribution function instead of the specification of the conformation of each individual molecule (in real systems there are too many molecules). Thus, the system can be described by a function $\Psi(\mathbf{x}, t, \mathbf{q}_1, \dots, \mathbf{q}_d)$ given the fraction of molecules that at position $\mathbf{x}$ and time $t$ have a conformation expressed by the vectors $\mathbf{q}_1, \dots, \mathbf{q}_d$. Now, the conservation of that distribution function results into an equation known as the Fokker-Planck equation whose solution is affected by the redoubtable *curse of dimensionality* because it is defined in a space of dimension $3+1+3 \cdot d$ (the physical space, the time and the $d$ conformation vectors). Even for a coarse description $d \approx 10$, the model is defined in a space involving 34 dimensions!

Other branches in science suffer this illness as for example the chemistry when the number of molecules of the reacting species is too low. In that case, the continuum approach based on the introduction of concentrations is no more a valid alternative. This situation is usually found in the modeling of genetic processes (e.g. expression of genes). Now, the state of the system is defined by the number of individuals of each one of the $d$ coexisting species $\Psi(\mathbf{x}, t, \#S_1, \dots, \#S_d)$ and again the balance equation governing the evolution of such system, the so-called chemical master equation, is defined in a multi-dimensional space that makes impossible its treatment by using the standard mesh based discretization techniques.

There are many other models defined in multidimensional spaces but for our discussion below these examples suffice.

In that concerns the issue related to the curse of dimensionality different techniques have been proposed for circumventing it, being Monte Carlo simulations the most widely used. The main drawback of stochastic simulations lies in the inevitable statistical noise affecting the solution, whose impact is important in particular when other than the moments of the distribution functions are computed [39]. Other possibility lies in the use of sparse grids based discretizations [16], within the deterministic framework, but they fail when the dimension of the space increases beyond a certain value (of about 20 dimensions as argued in [1]). In our knowledge there are few precedents of deterministic techniques able to circumvent efficiently the curse of dimensionality in highly multidimensional spaces. An appealing alternative that the author introduced some years ago in [6] for steady simulations and in [8] for transient solutions consists of 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) \times \dots \times F_i^d(x_d) \qquad (247)$$
*Remark.* In this expression the coordinates $x_i$ denote any coordinate, scalar or vectorial, involving the physical space, the time or any other conformation coordinate. □

Thus, if $M$ nodes are used to discretize each coordinate, the total number of unknowns involved in the solution is $N \times d \times M$ instead of the $M^d$ degrees of freedom involved in mesh based discretizations. We must recall that these functions are not "a priori" known, they are computed "on the fly" by introducing the approximation separated representation (247) into the model weak form and then solving the resulting non-linear problem. The interested reader can refer to [21] for a detailed description of the numerical and algorithmic aspects. The construction of such approximation is called Proper Generalized Decomposition because this decomposition is not orthogonal but for many usual transient models the number of terms in the finite sum is very close to the optimal decomposition obtained by applying the Proper Orthogonal Decomposition -POD- (or the Singular Value Decomposition -SVD-) on the model solution. As it can be noticed in the expression of the approximation separated representation (247) the complexity scales linearly with the dimension of the space in which the model is defined, instead of the exponential growing characteristic of mesh based discretization strategies. In general, for many models, the number of terms $N$ in the finite sum is quite reduced (few tens) and in all cases the approximation converges towards the solution associated with a fully tensorial product of the approximation bases considered in each space $x_i$. Thus, we can conclude about the generality of the separated representation, but its optimality depends on the solution features.

When we consider multidimensional models involving the physical space, the time, and a number of "exotic" extra-coordinates, the verdict concerning the CPU time savings when using the PGD is implacable: the PGD allows solving models never until now solved, suffering the so-called curse of dimensionality, and that were qualified many times as irresolvable! The PGD allows solving them, in some minutes, using a simple laptop.

Separated representations were applied for solving the multidimensional Fokker-Planck equation describing complex fluids in the kinetic theory framework. In [6], the authors addressed the solution of the linear Fokker-Planck equation describing multi-bead-spring molecular models in the steady state and in homogeneous flows. In that case the distribution function describing the molecular conformation only depends on the molecule configurational coordinates. The solution procedure was extended to non-linear kinetic theory formulations of more complex molecular descriptions in [48]. The transient solution was addressed in [8] in which the time was added as an extra-coordinate. Transient solutions making use of reduced bases in the context of an adaptive proper orthogonal decomposition [53] were considered in the case of low dimensional configuration spaces: the FENE model was addressed in [7] and liquid crystalline polymers in [10]. In [44], the authors considered multi-bead-spring models but considered spectral approximation for representing all the functions involved in the finite sums decomposition. A deeper analysis of non-linear and transient models was considered in [11]. Complex fluid models were coupled with complex flows in [52] and [49] opening very encouraging perspectives. In [20] the PGD was applied for solving the stochastic equation within the Brownian Configuration Field framework. Parametric models were addressed in [12]. The interested reader can refer to [21] [23] and the references therein for an exhaustive overview of PGD in computational rheology.

Multidimensional models encountered in the finer descriptions of matter (ranging from quantum chemistry to statistical mechanics descriptions) were revisited in [9]. The
multidimensional chemical master equation was efficiently solved in [22]. The Langer's equation governing phase transitions was considered in [41].

The application of separated representations for solving the highly multidimensional models describing concentrated suspensions of clusters constitutes a work in progress.

## References

1. Y. Achdou, O. Pironneau, Computational methods for option pricing, Siam Frontiers in Applied Mathematics, 2005.

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

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

4. S. Advani (ed.), Flow and Rheology in Polymer Composites Manufacturing, Elsevier, 1994.

5. A. Ammar, F. Chinesta, A particle strategy for solving the Fokker-Planck equation governing the fibre orientation distribution in steady recirculating flows involving short fibre suspensions, in *Lectures Notes on Computational Science and Engineering*, Springer, **43**, 1-16, 2005.

6. 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**, 153-176, 2006.

7. A. Ammar, D. Ryckelynck, F. Chinesta, R. Keunings, On the reduction of kinetic theory models related to finitely extensible dumbbells, *Journal of Non-Newtonian Fluid Mechanics*, **134**, 136-147, 2006.

8. 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**, 98-121, 2007.

9. A. Ammar, F. Chinesta, P. Joyot, The nanometric and micrometric scales of the structure and mechanics of materials revisited: An introduction to the challenges of fully deterministic numerical descriptions, *International Journal for Multiscale Computational Engineering*, 6/3, 191-213, 2008.

10. A. Ammar, E. Pruliere, F. Chinesta, M. Laso, Reduced numerical modeling of flows involving liquid-crystalline polymeres. *Journal of Non-Newtonian Fluid Mechanics*, **160**, 140-156, 2009.

11. A. Ammar, M. Normandin, F. Daim, D. Gonzalez, E. Cueto, F. Chinesta, Non-incremental strategies based on separated representations: Applications in computational rheology, *Communications in Mathematical Sciences*, **8/3**, 671-695, 2010.

12. A. Ammar, M. Normandin, F. Chinesta, Solving parametric complex fluids models in rheometric flows, *Journal of Non-Newtonian Fluid Mechanics*, **165**, 1588-1601, 2010.

13. G.K. Batchelor, The stress system in a suspension of force-free particles, J. Fluid Mech. **41**, 545-570, 1970.

14. N. Bellomo, Modeling complex living systems, Birkhauser, 2008.

15. R.B. Bird, C.F. Crutiss, R.C. Armstrong, O. Hassager, Dynamic of polymeric liquid, Volume 2: Kinetic Theory, John Wiley and Sons, 1987.

16. H.J. Bungartz, M. Griebel, Sparse grids, Acta Numerica, **13**, 1-123, 2004.

17. E. Cancès, M. Defranceschi, W. Kutzelnigg, C. Le Bris, Y. Maday, Computational Quantum Chemistry: a primer, *Handbook of Numerical Analysis*, Vol. X, Elsevier, 3-270, 2003.

18. C. Chaubal, A. Srinivasan, O. Egecioglu, G. Leal, Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems, J. Non-Newtonian Fluid Mech., **70**, 125-154, 1997.

19. 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**, 97-125, 2003.

20. F. Chinesta, A. Ammar, A. Falco, M. Laso. On the reduction of stochastic kinetic theory models of complex fluids, *Modeling and Simulation in Materials Science and Engineering*, **15**, 639-652, 2007.
21. F. Chinesta, A. Ammar, E. Cueto, Recent advances and new challenges in the use of the Proper Generalized Decomposition for solving multidimensional models, Archives of Computational Methods in Engineering, **17/4**, 327-350, 2010.

22. F. Chinesta, A. Ammar, E. Cueto, On the use of Proper Generalized Decompositions for solving the multidimensional Chemical Master Equation, European Journal of Computational Mechanics, **19**, 53-64, 2010.

23. F. Chinesta, A. Ammar, A. Leygue, R. Keunings, An overview of the Proper Generalized Decomposition with applications in computational rheology, Journal of Non Newtonian Fluid Mechanics, **166**, 578-592, 2011.

24. C. Cruz, L. Illoul, F. Chinesta, G. Regnier, Effects of a bent structure on the linear viscoelastic response of Carbon Nanotube diluted suspensions, Rheologica Acta, **49**, 1141-1155, 2010.

25. C. Cruz, F. Chinesta, G. Regnier, Review on the Brownian dynamics simulation of bead-rod-spring models encountered in computational rheology, Archives of Computational Methods in Engineering, **19/2**, 227-259, 2012.

26. E. Cueto, R. Monge, F. Chinesta, A. Poitou, I. Alfaro, M. Mackley, Rheological modeling and forming process simulation of CNT nanocomposites, International Journal of Material Forming, **3/2**, 1327-1338, 2010.

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

28. P. Dumont, S. Le Corre, L. Orgeas, D. Favier, A numerical analysis of the evolution of bundle orientation in concentrated fibre-bundle suspensions, Journal of Non-Newtonian Fluid Mechanics, **160**, 76-92, 2009.

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

30. J. Ferec, G. Ausias, M.C. Heuzey, P. Carreau, Modeling fiber interactions in semiconcentrated fiber suspensions, Journal of Rheology, **53/1**, 49-72, 2009.

31. F. Folgar, Ch. Tucker, Orientation behavior of fibers in concentrated suspensions, J. Reinf. Plast. Comp., **3**, 98-119, 1984.

32. P. Grassia, J. Hinch, L.C. Nitsche, Computer simulations of brownian motion of complex systems, J Fluid Mech., **282**, 373-403, 1995.

33. P. Grassia, J. Hinch, Computer simulations of polymer chain relaxation via brownian motion, J Fluid Mech., **308**, 255-288, 1996.

34. G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech. **13**, 33-62, 1962.

35. J. Hinch, G. Leal, The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles, J. Fluid Mech., **52**, 683-712, 1972.

36. J. Hinch, G. Leal, Constitutive equations in suspension mechanics. Part I, J. Fluid Mech., **71**, 481-495, 1975.

37. J. Hinch, G. Leal, Constitutive equations in suspension mechanics. Part II, J. Fluid Mech., **76**, 187-208, 1976.

38. G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proc. R. Soc. London, **A102**, 161-179, 1922.

39. R. Keunings, Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory, Rheology Reviews, D.M. Binding and K. Walters (Edts.), British Society of Rheology, 67-98, 2004.

40. M. Kroger, A. Ammar, F. Chinesta, Consistent closure schemes for statistical models of anisotropic fluids, Journal of Non-Newtonian Fluid Mechanics, **149**, 40-55, 2008.

41. H. Lamari, A. Ammar, A. Leygue, F. Chinesta, On the solution of the multidimensional Langer's equation by using the Proper Generalized Decomposition Method for modeling phase transitions, Modelling and Simulation in Materials Science and Engineering, **20**, 015007 (2012).

42. R.B. Laughlin, The theory of everything, Proceeding of the U.S.A. National Academy of Sciences, 2000.

43. C. Le Bris (ed.), Handbook of Numerical Analysis, Vol. X: Computational Chemistry, Elsevier, 2003.

44. G.M. Leonenko, T.N. Phillips, On the solution of the Fokker-Planck equation using a high-order reduced basis approximation, Computer Methods in Applied Mechanics and Engineering **199/1-4**, 158-168, 2009.

45. A. Ma, F. Chinesta, M. Mackley, A. Ammar, The rheological modelling of Carbon Nanotube (CNT) suspensions in steady shear flows, International Journal of Material Forming, **2**, 83-88, 2008.
46. A. Ma, F. Chinesta, A. Ammar, M. Mackley, Rheological modelling of Carbon Nanotube aggregate suspensions, Journal of Rheology, **52/6**, 1311-1330, 2008.

47. A. Ma, F. Chinesta, M. Mackley, The rheology and modelling of chemically treated Carbon Nanotube suspensions, Journal of Rheology, **53/3**, 547-573, 2009.

48. B. Mokdad, E. Pruliere, A. Ammar, F. Chinesta, On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach, Applied Rheology, **17/2**, 26494, 1-14, 2007.

49. B. Mokdad, A. Ammar, M. Normandin, F. Chinesta, J.R. Clermont, A fully deterministic micro-macro simulation of complex flows involving reversible network fluid models, Mathematics and Computer in Simulation, **80**, 1936-1961, 2010.

50. H.C. Ottinger, Stochastic processes in polymeric fluids. Springer, Berlin Heidelberg, 1996.

51. C. Petrie, The rheology of fibre suspensions, J. Non-Newtonian Fluid Mech., **87**, 369-402, 1999.

52. E. Pruliere, A. Ammar, N. El Kissi, F. Chinesta, Recirculating flows involving short fiber suspensions: Numerical difficulties and efficient advanced micro-macro solvers, Archives of Computational Methods in Engineering, State of the Art Reviews, **16**, 1-30, 2009.

53. D. Ryckelynck, F. Chinesta, E. Cueto, A. Ammar, On the a priori model reduction: Overview and recent developments, Archives of Computational Methods in Engineering, State of the Art Reviews, **13/1**, 91-128, 2006.

54. Ch. Tucker, Flow regimes for fiber suspensions in narrow gaps, J. Non-Newtonian Fluid Mech., **39**, 239-268, 1991.