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

HAL
open science

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

Marta Perez, Adrien Scheuer, Emmanuelle Abisset-Chavanne, Francisco Chinesta, Roland Keunings

► To cite this version:

Marta Perez, Adrien Scheuer, Emmanuelle Abisset-Chavanne, Francisco Chinesta, Roland Keunings. A multi-scale description of orientation in simple shear flows of confined rod suspensions. Journal of Non-Newtonian Fluid Mechanics, 2016, 233, pp.61 - 74. 10.1016/j.jnnfm.2016.01.011. hal-01403860

HAL Id: hal-01403860

https://hal.science/hal-01403860

Submitted on 28 Nov 2016

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
# A multi-scale description of orientation in simple shear flows of confined rod suspensions

M. Perez<sup>a</sup>, A. Scheuer<sup>a,b</sup>, E. Abisset-Chavanne<sup>a</sup>, F. Chinesta<sup>a</sup>, R. Keunings<sup>b</sup>

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

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

The multi-scale description of dilute or semi-dilute suspensions involving rods has been successfully accomplished and applied in many scenarios of industrial interest. Many processes involve, however, the flow of rod suspensions in very narrow gaps whose thickness is much smaller than the rod length. In these conditions, the evolution of rod orientation is expected to be affected by confinement effects. In the present work, we propose a multi-scale description of rod orientation in confined conditions and simple shear flows.

## 1. Introduction

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

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

Mesoscopic kinetic theory models result from coarsening microscopic descriptions. In kinetic theory models, the individuality of the particles is lost in favour of a statistical description that substitutes the entities by a series of conformation coordinates [6,12]. For example, when considering a suspension of rods, the

mesoscopic description consists in giving the fraction of rods that at position **x** and time *t* are oriented along direction **p**. This information is contained in a probability distribution function – pdf – whose evolution is governed by a Fokker-Planck equation. Fokker-Planck equations being multidimensional (in time, physical space and conformation space), standard mesh-based discretization techniques fail when addressing their numerical solution. This issue is known as the "curse of dimensionality". The direct solution of Fokker-Planck equations has been made feasible in many cases with the introduction of the Proper Generalized Decomposition (PGD) approach [8,9].

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

In the case of dilute suspensions of short fibres in a Newtonian fluid, the three scales have been extensively considered without major difficulties to model the associated systems. Challenges appear, however, as soon as the concentration increases. In the semi-dilute and semi-concentrated regimes, fibre-fibre interactions occur, but in general they can be accurately modelled by introducing a randomizing diffusion term [15]. There is a wide literature on the modelling of dilute and semi-dilute suspensions, e.g. [4,17–20,27]. Available models describe quite well the experimental observations.

When fibre orientation predictions are compared with experimental results in injection processes involving concentrated suspensions of short fibres, a noticeable delay in the orientation kinematics is observed. Some ad hoc modelling approaches were
proposed to delay the orientation kinematics predicted with the
Folgar and Tucker model [15]. In these models, a fluid-particle slid-
ing mechanism is introduced or different rate equations for the
eigenvalues and eigenvectors of the orientation tensor are used
[28,31,32]. A more physically-based approach was proposed by
Ferec and co-authors in [14], where the interaction mechanisms
were taken into account within a multi-scale framework.

All these studies concerned unconfined flows, despite the fact
that processes of industrial interest often involve narrow gaps
where fibres may have length greater than the gap and wall ef-
fects cannot be ignored. The orientation delay was indeed observed
in such confined flows (e.g. [24]), wherein the orientation process
could probably differ from the one predicted by the standard Jeff-
ery equation and the mesoscopic and macroscopic models derived
from it.

In the present paper, we propose a multi-scale description of rod orientation in confined conditions and analyse the impact of initial conditions on rod kinematics under confinement constraints. It is well known that the orientation kinematics of fibres located near the wall deviate from Jeffery's predictions, even when the fibres never enter in contact with the wall [16,21,25,30]. The model proposed here does not take into account these hydrodynamic confinement effects, but rather it considers physical contact between rods and gap walls through the introduction of the contact force ensuring wall impenetrability. In order to consider the rod kinematics perturbation in absence of physical contact with the wall, we should consider an extra hydrodynamic force acting on the dumbbell beads describing the rod (see below); this force would depend on the distance to the wall and the approaching velocity. The model proposed in this work aims at capturing the first-order effects of confinement. The effect of the wall proximity on the fiber kinematics constitutes a second-order effect, to be addressed in future works.

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

* if **a** and **b** are first-order tensors, then the single contraction **"·"** reads (a · b) = a_j b_j;

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

* if **a** and **b** are respectively second and first-order tensors, then the single contraction **"·"** reads (a · b)_j = a_{jk} b_k;

* if **a** and **b** are respectively third and first-order tensors, then the single contraction **"·"** reads (a · b)_{jk} = a_{jkm} b_m;

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

* if **a** and **b** are respectively second and fourth-order tensors, then the double contraction **":**" reads (a : b)_{jk} = a_{ml} b_{mljk}.

## 2. Modelling confined suspensions of rods

We consider a dilute suspension of rigid, non-Brownian, high-
aspect-ratio fibres suspended in a Newtonian fluid of viscosity η.
The fibres are modelled as rigid rods of length 2L. The rod orienta-
tion is given by the unit vector **p** located at the rod centre of grav-
ity G and aligned along the rod axis. Inertial effects are neglected
in the sequel. We assume that the presence and orientation of the
rods do not affect the flow kinematics that is defined by the ve-
locity field **v**^T = (γz, 0, 0) describing a simple shear flow. Elonga-
tion is not considered in this work because it tends to reduce and
even suppress all confinement effects when its intensity with re-
spect to shear is large enough. In the case of pure elongation, the
rods align monotonically in the extension direction and can never
enter in contact with the gap walls. The flow occurs in a narrow
gap Ω × [-H, H], with x^T = (x, y) ∈ Ω ⊂ R^2 assumed large enough
and z ∈ [-H, H] with H < L for ensuring confinement conditions.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Hydrodynamic and contact forces acting on a confined rod immersed in a simple shear flow.</figcaption></figure>

We consider first an individual rod whose centre of gravity is lo-
cated in the mid plane z_G = 0. Thus, both rod extremities enter in
contact or lose contact with the gap walls simultaneously. With-
out loss of generality, we consider that the unit vector **p** related
to bead located at **p**L points towards the upper wall, and conse-
quently **-p** points towards the opposite wall.

As just indicated, this work only considers confined flow ex-
hibiting constant shear rate throughout the gap thickness in order
to ensure a constant shear rate along the rod length. This flow
is conceptually of interest because it allows one to identify the
differences between the standard (unconfined) Jeffery model and
the one that results when confinement effects take place. It re-
mains, however, quite far from the applicative processing condi-
tions that in general involve more complex flows, in many cases of
Poiseuille type. In processing conditions, almost-parabolic velocity
profiles throughout the gap thickness are usually encountered, im-
plying that the velocity gradient is no longer constant along the
length of the rod when its extremities approach both walls. In
these circumstances, one should consider second-order kinemat-
ics in the derivation of the rod rotary velocity, as done in [2,7] for
unconfined flows. In that applicative perspective, one should also
take into account non-Newtonian rheology of the suspending fluid,
which adds a major difficulty as discussed in [7]. All these issues
will be presented and discussed in ongoing publications. In the
present work, we concentrate on simple shear flows, Newtonian
suspending fluids, and we assume incompressible and isothermal
conditions.

2.1. Microscopic description of confined kinematics of an individual rod

It is well known that the kinematics of a rod with infinite as-
pect ratio immersed in an unconfined simple shear flow with ve-
locity gradient ∇**v** is given by Jeffery's equation [22]

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

Jeffery's equation predicts full alignment in the flow direc-
tion. For finite aspect ratio ellipsoidal particles, periodic trajecto-
ries known as Jeffery's orbits are predicted instead. Brownian ef-
fects also avoid full alignment and are generally introduced at the
rod population level.

When the domain thickness is smaller than the rod length (H < L), some Jeffery's trajectories are forbidden, i.e. those trajectories involving Lp(t) · n > H, where n = (0, 0, 1)^T is the unit vector defining the thickness direction. In that case, the rod kinematics are defined by the standard Jeffery model (1) while p · n < H/L and are perturbed as soon as the rod reaches the upper and lower walls. In order to determine the perturbed rod kinematics, we represent the rod as a dumbbell [1] with hydrodynamic and contact forces acting on the dumbbell beads as illustrated in Fig. 1. We assume that hydrodynamic forces applied on each bead F^H depend
<figure><img src="image_3.png" /><figcaption>Fig. 2. Orientation evolution for H/L = 0.8. Confined (red curve) versus unconfined (broken blue curve) trajectories.</figcaption></figure>

on the difference of velocities between the fluid and the bead, the first one given by $\mathbf{v}_0 + \nabla\mathbf{v} \cdot \mathbf{pL}$ and the second one by $\mathbf{v}_G + \dot{\mathbf{p}}L$. Thus, the force $\mathbf{F}^H(\mathbf{pL})$ is given by

$$ \mathbf{F}^H(\mathbf{pL}) = \xi (\mathbf{v}_0 + \nabla\mathbf{v} \cdot \mathbf{pL} - \mathbf{v}_G - \dot{\mathbf{p}}L), \quad (2) $$

where $\xi$ is the friction coefficient, $\mathbf{v}_0$ the fluid velocity at the rod centre of gravity, and $\mathbf{v}_G$ the velocity of the centre of gravity. The contact force is assumed to act in the direction perpendicular to the wall (we thus ignore possible friction at the wall):

$$ \mathbf{F}^C(\mathbf{pL}) = \mu \mathbf{n}, \quad (3) $$

with $\mathbf{F}^C(\mathbf{pL}) = -\mathbf{F}^C(-\mathbf{pL})$. The contact force intensity $\mu$ is unknown at this stage and will be determined such as to prevent the beads from penetrating the walls.

Thus, the total force acting on the bead located at $\mathbf{pL}$ is $\mathbf{F}(\mathbf{pL}) = \mathbf{F}^H(\mathbf{pL}) + \mathbf{F}^C(\mathbf{pL})$. As inertia is neglected, the resultant force acting on the dumbbell vanishes. Since $\mathbf{F}^C(\mathbf{pL}) = -\mathbf{F}^C(-\mathbf{pL})$, the force balance yields $\mathbf{F}^H(\mathbf{pL}) = -\mathbf{F}^H(-\mathbf{pL})$ and thus

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

i.e. the velocity of the rod centre of gravity coincides with the velocity of the fluid at that position. This ensures that both beads are simultaneously in contact with the walls.

The resulting torque must also vanish, that is $\mathbf{F} = \mathbf{F}^H + \mathbf{F}^C = \lambda \mathbf{p}$, with $\lambda \in \mathbb{R}$. Since $\mathbf{v}_0 = \mathbf{v}_G$, we have

$$ \xi (\nabla\mathbf{v} \cdot \mathbf{pL} - \dot{\mathbf{p}}L) + \mu \mathbf{n} = \lambda \mathbf{p}. \quad (5) $$

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

$$ \lambda = \xi L (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) + \mu p_z, \quad (6) $$

with $p_z = \mathbf{p} \cdot \mathbf{n}$. Injecting this expression for $\lambda$ into Eq. (5) yields

$$ \xi L (\nabla \mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}) + \mu \mathbf{n} = \xi L (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p} + \mu p_z \mathbf{p}, \quad (7) $$

or

$$ (\nabla \mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}) + \frac{\mu}{\xi L} \mathbf{n} = (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p} + \frac{\mu}{\xi L} p_z \mathbf{p}. \quad (8) $$

The rotary velocity $\dot{\mathbf{p}}$ is thus given by

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p} + \frac{\mu}{\xi L} (\mathbf{n} - p_z \mathbf{p}) = \dot{\mathbf{p}}^J + \dot{\mathbf{p}}^C. \quad (9) $$

Here, $\dot{\mathbf{p}}^J$ denotes the Jeffery rotary velocity component given by Eq. (1), while the confined component is defined as $\dot{\mathbf{p}}^C = \frac{\mu}{\xi L} (\mathbf{n} - p_z \mathbf{p})$.

We now need to determine the contact force intensity $\mu$. In order to obtain its value we must consider that the contact force appears in order to avoid that the bead leaves the flow domain, that is $\dot{\mathbf{p}} \cdot \mathbf{n} \le 0$, with $\mu \ne 0$ if $p_zL = H$ and $\dot{\mathbf{p}}^J \cdot \mathbf{n} > 0$. The contact force $\mu \mathbf{n}$ must ensure that the resulting velocity is tangent to the wall. Because dynamical effects are neglected, the bead cannot rebound. It thus suffices to enforce the condition $\dot{\mathbf{p}} \cdot \mathbf{n} = 0$. Multiplying Eq. (9) by $\mathbf{n}$ yields

$$ 0 = \dot{\mathbf{p}}^J \cdot \mathbf{n} + \frac{\mu}{\xi L} (1 - p_z^2) = [\dot{\mathbf{p}}^J]_z + \frac{\mu}{\xi L} (1 - p_z^2), \quad (10) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Orientation evolution for H/L = 0.8 (top-left), H/L = 0.6 (top-right), H/L = 0.4 (bottom-left) and H/L = 0.2 (bottom-right).</figcaption></figure>

or

$$ \mu = - \frac{\xi L}{1 - p_z^2} [\dot{\mathbf{p}}^J]_z. \qquad (11) $$

Eq. (9) thus reduces to

$$ \dot{\mathbf{p}} = \dot{\mathbf{p}}^J - \frac{1}{1 - p_z^2} [\dot{\mathbf{p}}^J]_z (\mathbf{n} - p_z \mathbf{p}) = \dot{\mathbf{p}}^J + \dot{\mathbf{p}}^C, \qquad (12) $$

with $\dot{\mathbf{p}}^J$ given by Eq. (1), $p_z = \mathbf{p} \cdot \mathbf{n}$ and $[\dot{\mathbf{p}}^J]_z = \dot{\mathbf{p}}^J \cdot \mathbf{n}$.

The first term in Eq. (12) corresponds to the standard, uncon- fined Jeffery kinematics, while the second term avoids that rod beads leave the flow domain. While preparing this paper, it came to our attention that Eq. (12) had been derived independently in [26] using similar arguments but a different procedure.

Eq. (9) can be rewritten as

$$ \dot{\mathbf{p}} - \mathbf{\Omega} \cdot \mathbf{p} = \mathbf{D} \cdot \mathbf{p} - (\mathbf{D} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p} + \frac{\mu}{\xi L} (\mathbf{n} - p_z \mathbf{p}), \qquad (13) $$

from which we can conclude that the rate of strain due to the flow kinematics and the reaction forces induce an effective rotation $\dot{\mathbf{p}} - \mathbf{\Omega} \cdot \mathbf{p}$, thus ensuring the model objectivity.

So, to sum up, the kinematics of a rigid rod of length 2L in a simple shear flow occurring in a narrow gap of thickness 2H are

given by the following confined Jeffery model:

$$ \left\{ \begin{array}{ll} \dot{\mathbf{p}} = \dot{\mathbf{p}}^J & \text{if } p_z L < H \\ \dot{\mathbf{p}} = \dot{\mathbf{p}}^J & \text{if } p_z L = H \text{ & } \dot{\mathbf{p}}^J \cdot \mathbf{n} \le 0, \\ \dot{\mathbf{p}} = \dot{\mathbf{p}}^J + \dot{\mathbf{p}}^C & \text{if } p_z L = H \text{ & } \dot{\mathbf{p}}^J \cdot \mathbf{n} > 0 \end{array} \right. \qquad (14) $$

where the rods are in contact with the walls if $p_z = H/L$. It is im- portant to notice that expression (12) only applies if the trial Jeff- fery velocity is such that $\dot{\mathbf{p}}^J \cdot \mathbf{n} > 0$. When $p_z L < H$, $\dot{\mathbf{p}}$ reduces to the standard Jeffery contribution $\dot{\mathbf{p}}^J$. When $p_z L = H$ and $\dot{\mathbf{p}}^J \cdot \mathbf{n} \le 0$, the rotary velocity also reduces to the Jeffery contribution. The contact being unilateral, a rod can detach from the wall as soon as its ve- locity induces the detachment. These conditions, known as Kuhn and Tucker or Signorini conditions, are similar to those encoun- tered when describing unilateral contact or elastoplasticity. When projecting Eq. (12) in the thickness direction, we have $\dot{\mathbf{p}}^J \cdot \mathbf{n} = 0$, thus ensuring that the rod beads do not leave the flow domain.

In the numerical tests discussed below, we will occasionally as- sume fully-confined conditions, meaning that fibres are assumed to be and always remain in contact with the walls, while their orien- tation evolution is governed by Eq. (12). In such case, the contact is bilateral and fibres are prevented to detach from the wall even when $\dot{\mathbf{p}}^J \cdot \mathbf{n} \le 0$. A reaction force must appear to maintain contact
<figure><img src="image_5.png" /><figcaption>Fig. 4. Confined Jeffery solution $a^{\text{disc}}(t)$ versus fully-confined macroscopic orientation tensor $a^c(t)$ in the case of rods aligned in the same direction. It can be noticed that the planar components of the orientation tensor evolve smoothly whereas the out-of-plane component shows a localized evolution when, after a period of time during which the fibre beads slide at the wall, they detach from it.</figcaption></figure>

and avoid detachment of the rod bead from the wall. In the sequel, this particular model will be referred to as the fully-confined Jeffery model. Although not entirely physical, it will be useful for validation purposes.

## 2.2. Mesoscopic description of a population of rods

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

### 2.2.1. Discrete description

The discrete approach consists in computing the orientation of each individual rod belonging to a large discrete ensemble of $N$ rods. Thus, the population is described from the individuals composing it, whose conformation is given by vectors $\mathbf{p}_i$, $i = 1, ..., N$,

each governed by Eq. (14). The main drawback of this approach lies in the necessity of tracking the evolution of each rod by solving the corresponding equation, and even if conceptually there is no major difficulty, the computing cost could be excessive in most practical applications.

### 2.2.2. Continuous description

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

The pdf satisfies the normalization condition:

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

where $S$ is the rod conformation space, i.e. the surface of the unit sphere in the unconfined case.

Conservation of probability leads to the so-called Fokker-Planck equation

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

where the rod rotary velocity $\dot{\mathbf{p}}$ is given by the confined Jeffery's equation (14) and $\dot{\mathbf{x}} = \mathbf{v}(\mathbf{x}, t)$.

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

## 2.3. Macroscopic description

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

<figure><img src="image_6.png" /><figcaption>Fig. 5. Fully-confined Jeffery solution $a^{\text{disc,fc}}$ (fibre detachment from the wall is prevented) versus fully-confined macroscopic description $a^c(t)$ in the case of rods aligned in the same direction.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Isotropic initial orientation distribution for $H/L = 0.2$ and $N = 4000$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Evolution of the confined (discontinuous line) and unconfined (continuous line) components of the second-order orientation tensor for $H/L = 0.2$.</figcaption></figure>

### 2.3.1. Discrete approach

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

$$ \mathbf{a}^{\text{disc}}(t) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{p}_i(t) \otimes \mathbf{p}_i(t). \qquad (17) $$

Here, the superscript disc refers to the discrete approach.

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

### 2.3.2. Continuous approach: the two limiting cases

In the continuous approach, the orientation distribution function is substituted by its moments for describing the microstructure [3] at the macroscale, and an evolution equation for these moments is derived from the Fokker-Planck equation. Usually, macroscopic descriptions of rod suspensions are based on the use of the first two non-vanishing moments, i.e. the second and fourth-order moments, **a** and **A**, respectively defined by

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

and

$$ \mathbf{A} = \int_S \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p}. \qquad (19) $$
<figure><img src="image_9.png" /><figcaption>Fig. 8. Initial in-plane isotropic orientation distribution related to the fully-confined simulation on z = H, H/L = 0.2, implying p_z = 0.2 for all fibres.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Confined Jeffery's orientation tensor $\mathbf{a}^{\text{disc}}(t)$ versus fully-confined macroscopic description $\mathbf{a}^c(t)$ for a population of rods with the fully-confined initial isotropic distribution depicted in Fig. 8.</figcaption></figure>

Odd moments vanish in view of the symmetry of the pdf: $\psi(\mathbf{x}, t, \mathbf{p}) = \psi(\mathbf{x}, t, -\mathbf{p})$.

We consider two extreme situations, namely the unconfined case which yields the standard macroscopic description based on Jeffery's equation, and the fully-confined case wherein all rods are in contact with the gap walls. The latter is certainly relevant in conditions of intense confinement, namely for $H \ll L$. In this case, we can further assume that rod orientations are distributed in the allowed region of the unit sphere but under a lubrication kinematic constraint.

These three modelling frameworks are described below.

* The unconfined case: For unconfined rods, the orientation is defined on the surface of the unit sphere $S$, and the associated Fokker-Planck equation is readily exploited to derive the classical evolution equation for the second-order orientation tensor associated to the standard Jeffery model [5]:

$$ \dot{\mathbf{a}}^J = \nabla \mathbf{v} \cdot \mathbf{a}^J + \mathbf{a}^J \cdot (\nabla \mathbf{v})^T - 2 \mathbf{A}^J : \nabla \mathbf{v}. \quad (20) $$

Here, the superscript $J$ indicates that we consider the orientation tensors associated with the standard Jeffery model for unconfined systems.

* The fully-confined case: In the other limiting case, we assume that all rods in the suspension are in contact with the walls and that the contact is bilateral, i.e. rods can orient but always remain in contact with the walls. In this situation, we can derive an evolution equation for what we call the fully-confined orientation tensor $\mathbf{a}^c$ defined as

$$ \mathbf{a}^c = \int_C \mathbf{p} \otimes \mathbf{p} \, \psi(\mathbf{p}) \, d\mathbf{p}, \quad (21) $$

with $C = \{\mathbf{p}, p_z = \pm H/L\}$.

Indeed, we have

$$ \dot{\mathbf{a}}^c = \int_C \mathbf{p} \otimes \mathbf{p} \, \dot{\psi}(\mathbf{p}) \, d\mathbf{p}, \quad (22) $$

where

$$ \dot{\psi} = -\nabla_p \cdot (\dot{\mathbf{p}} \, \psi(\mathbf{p})). \quad (23) $$

Integration by parts with respect to coordinates $\mathbf{p}$ yields

$$ \dot{\mathbf{a}}^c = \int_C (\dot{\mathbf{p}} \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}) \, \psi(\mathbf{p}) \, d\mathbf{p}. \quad (24) $$

In the confined case

$$ \mathbf{p} = \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix}. \quad (25) $$
<figure><img src="image_11.png" /><figcaption>Fig. 10. Fully-confined Jeffery's orientation tensor versus $a^C(t)$ for a population of rods with the fully-confined initial isotropic distribution depicted in Fig. 8.</figcaption></figure>

with $p_z = \pm H/L$, from which Eq. (21) gives

$$ \mathbf{a}^C = \int_C \begin{pmatrix} \mathbf{q} \otimes \mathbf{q} & p_z \mathbf{q} \\ p_z \mathbf{q}^T & p_z^2 \end{pmatrix} \psi(\mathbf{p}) \, d\mathbf{p}. \qquad (26) $$

Considering now Eq. (24), with

$$ \dot{\mathbf{p}} \otimes \mathbf{p} = \begin{pmatrix} \dot{\mathbf{q}} \\ 0 \end{pmatrix} \otimes \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix}, \qquad (27) $$

and

$$ \mathbf{p} \otimes \dot{\mathbf{p}} = \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix} \otimes \begin{pmatrix} \dot{\mathbf{q}} \\ 0 \end{pmatrix}, \qquad (28) $$

we obtain

$$ \dot{\mathbf{a}}^C = \int_C \begin{pmatrix} (\dot{\mathbf{q}} \otimes \mathbf{q} + \mathbf{q} \otimes \dot{\mathbf{q}}) & p_z \dot{\mathbf{q}} \\ p_z \dot{\mathbf{q}}^T & 0 \end{pmatrix} \psi(\mathbf{p}) \, d\mathbf{p}. \qquad (29) $$

Now, we define the in-plane second-order orientation tensor $\mathbf{b}$ according to

$$ \mathbf{b} = \int_C \mathbf{q} \otimes \mathbf{q} \, \psi(\mathbf{p}) \, d\mathbf{p}. \qquad (30) $$

This new tensor does not have a unit trace since $||\mathbf{q}||^2 = 1 - \frac{H^2}{L^2}$. We also define a vector $\mathbf{c}$,

$$ \mathbf{c} = \int_C \mathbf{q} \, p_z \, \psi(\mathbf{p}) \, d\mathbf{p}, \qquad (31) $$

which does not vanish because the pdf is symmetric, i.e. $\psi(\mathbf{p}) = \psi(-\mathbf{p})$ and $\mathbf{q}$ and $p_z$ have opposite sign on each parallel $z = \pm H$ defining $\mathbf{c}$.

With these definitions and Eq. (26), we obtain the orientation tensor $\mathbf{a}^C$,

$$ \mathbf{a}^C = \begin{pmatrix} \mathbf{b} & \mathbf{c} \\ \mathbf{c}^T & \frac{H^2}{L^2} \end{pmatrix}, \qquad (32) $$

and its time derivative reads

$$ \dot{\mathbf{a}}^C = \begin{pmatrix} \dot{\mathbf{b}} & \dot{\mathbf{c}} \\ \mathbf{c}^T & 0 \end{pmatrix}, \qquad (33) $$

with $\dot{\mathbf{b}}$ and $\dot{\mathbf{c}}$ given by Eq. (29).

In order to obtain $\mathbf{b}$ and $\dot{\mathbf{c}}$, we need to derive the expression of $\dot{\mathbf{q}}$. For that purpose, we decompose the velocity gradient according to

$$ \nabla \mathbf{v} = \begin{pmatrix} \mathbf{G} & \mathbf{g} \\ \mathbf{j}^T & \mathcal{G} \end{pmatrix}, \qquad (34) $$

such that

$$ \nabla \mathbf{v} \cdot \mathbf{p} = \begin{pmatrix} \mathbf{G} \\ \mathbf{j}^T \end{pmatrix} \begin{pmatrix} \mathbf{g} \\ \mathcal{G} \end{pmatrix} \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix} = \begin{pmatrix} \mathbf{G} \cdot \mathbf{q} + \mathbf{g} p_z \\ \mathbf{j}^T \cdot \mathbf{q} + \mathcal{G} p_z \end{pmatrix}. \qquad (35) $$

As detailed in Appendix B, we obtain

$$ \begin{aligned} \dot{\mathbf{q}} = \mathbf{G} \cdot \mathbf{q} - \delta_1 (\mathbf{G} : (\mathbf{q} \otimes \mathbf{q})) \mathbf{q} &+ \delta_2 \mathbf{g} \\ & - \delta_3 (\mathbf{q}^T \cdot \mathbf{g}) \mathbf{q}, \end{aligned} \qquad (36) $$

with $\delta_1 = \frac{1}{1-p_z^2}$, $\delta_2 = p_z$ and $\delta_3 = \frac{p_z}{1-p_z^2} = \delta_1 \delta_2$.

Thus, the time derivative of tensor $\mathbf{b}$ is given by

$$ \begin{aligned} \dot{\mathbf{b}} = \mathbf{G} \cdot \mathbf{b} + \mathbf{b} \cdot \mathbf{G}^T - 2\delta_1 \mathbf{G} : \mathbf{B} \\ + (\mathbf{g} \otimes \mathbf{c} + \mathbf{c} \otimes \mathbf{g}) - 2\delta_3 \mathcal{B} \cdot \mathbf{g}. \end{aligned} \qquad (37) $$

with $\mathcal{B}$ and $\mathbf{B}$ being respectively the third and fourth-order in-plane orientation tensors related to the in-plane orientation $\mathbf{q}$. Both will be expressed from $\mathbf{c}$ and $\mathbf{b}$ through adequate closure relations. The simplest closure consists in assuming

$$ \begin{cases} \mathcal{B} = \frac{1}{p_z} \mathbf{b} \otimes \mathbf{c} \\ \mathbf{B} = \mathbf{b} \otimes \mathbf{b} \end{cases}. \qquad (38) $$

It is exact when the rods are fully aligned (i.e. when the pdf reduces to a Dirac delta distribution). In the general case, its validity must be checked carefully. The development of more accurate closure approximations following the rationale considered in [10,11,23,29] constitutes a work in progress that will be reported in future publications.

Finally, we obtain the time evolution of the first-order moment

$$ \dot{\mathbf{c}} = \int_C \dot{\mathbf{q}} \, p_z \, \psi(\mathbf{p}) \, d\mathbf{p}. \qquad (39) $$

by considering the expression of $\dot{\mathbf{q}}$ given by Eq. (36),

$$ \dot{\mathbf{c}} = \mathbf{G} \cdot \mathbf{c} - \delta_1 \delta_2 \mathbf{G} : \mathcal{B} + \delta_2^2 \mathbf{g} - \delta_3 \delta_2 \mathbf{b} \cdot \mathbf{g}. \qquad (40) $$

* The lubrication simplified model: Inspired by lubrication theory that successfully reduces the 3D flow equations in the case of thin gaps by neglecting out-of-plane velocities, in the case of intense confinement, i.e. $H/L < 0.3$, we could assume $\dot{\mathbf{p}}_z \approx 0$. Moreover, the in-plane components of $\dot{\mathbf{p}}^T$ and of its confined counterpart $\dot{\mathbf{p}}$ given by (12) are very close in such situation. Thus, one could ignore the unconfined motion of rods in the thin gap and consider that all rods are each fully confined by an imaginary wall located at the initial bead z-coordinates.

This is equivalent to computing the fully-confined solution $\mathbf{a}^C$ at each $z$, $\mathbf{a}^C(t; z)$, $z \in [0, H]$, and then considering the average

$$ \dot{\mathbf{a}}^C(t; H) = \int_0^H \omega(z) \mathbf{a}^C(t; z) dz, \qquad (41) $$

where $\omega(z)$ is the fraction of rods that initially have a bead located at $z$ (the other bead being at $-z$).

Thus, $\dot{\mathbf{a}}^C$ could be retained as a simplified orientation descriptor when confinement becomes dominant.

In what follows, $\dot{\mathbf{a}}^C$ is calculated by integrating $\mathbf{a}^C$ at different coordinates $z_i$ uniformly distributed in the gap $[0, H]$ and then averaging these solutions. For an isotropic initial distribution in the gap, we obtain the through-the-thickness average of $\mathbf{a}^C$.

Another possibility, not exploited in the present paper, consists in deriving the equation governing the time evolution of tensor $\dot{\mathbf{a}}^C$ by introducing Eqs. (37) and (40) into the time derivative of Eq. (41) and then performing integration. The derivation is detailed in Appendix C, where we show that, even though a closed evolution equation does not exist, suitable approximate expressions can be obtained.
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Comparison of confined $a^{disc}$ and unconfined $a^J_{disc}$ solutions computed by using a population of rods large enough to ensure their use as reference solutions, and the orientation tensor $\mathbf{a}^J(t)$ obtained by integrating Eq. (20) using the quadratic and hybrid closure relations, $a^{quad}$ and $a^{hybr}$ respectively, for different initial confined orientation distributions.</figcaption></figure>

### 3. Model predictions in simple shear flow

In this section, we discuss the predictions of the proposed models for the case of a simple shear flow with velocity $\mathbf{v}^T = (\dot{\gamma}z, 0, 0)$ and $\dot{\gamma} = 1$. Since we consider a unit shear rate, the time coordinate in most of the graphical representations that follow can be viewed as a shear strain coordinate.

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

#### 3.1. Evaluating the trajectory of a single rod

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

##### 3.1.1. Solution of the confined Jeffery equation

Here, we track the orientation of a rod initially unconfined, that is $p_zL < H$, but whose Jeffery trajectory implies at a certain instant

that $p_z^J(t)L > H$. It is important to note that in our model all initial orientations associated with Jeffery's trajectories never reaching the domain wall (or reaching it at the highest point where the trajectory becomes tangent to the wall) will follow a Jeffery trajectory without any perturbation. As indicated in the introduction, hydrodynamic effects appear when the rod approaches the gap walls. These second-order effects are ignored in the simulations that follow.

Fig. 2 shows the orientation trajectory followed by the rod. At the beginning, because the rod is unconfined, it follows the trajectory dictated by Jeffery's equation (1). The actual trajectory $\mathbf{p}(t)$ is depicted in red whereas the unconfined Jeffery motion $\mathbf{p}^J(t)$ is represented by the broken blue curve. Until reaching the walls, as expected, both trajectories superpose and consequently cannot be distinguished. As soon as one of the rod beads reaches the upper wall $p_zL = H$ (the other bead reaches simultaneously the lower wall, but by using symmetry arguments we only refer to the one touching the upper wall), with $\dot{\mathbf{p}}^J \cdot \mathbf{n} > 0$ for both trajectories, the actual trajectory and the one associated to the unconfined Jeffery motion bifurcate from each other. The Jeffery trajectory continues its unconfined motion with $p_z^JL > H$ whereas the confined rod continues to slide on the upper wall $p_zL = H$ until it detaches from it.

The orientations satisfying the contact condition $p_zL = H$ define the two thick black parallels depicted in Fig. 2. Thus, we notice that during a time interval the rod orients toward the flow direction while remaining in contact with the walls. In fact, during this period of time, the confined rod is crossing an infinity of forbidden Jeffery trajectories with $\dot{\mathbf{p}}^J \cdot \mathbf{n} > 0$. Finally, at a certain instant the rod reaches an unconfined Jeffery orbit that is tangent to the wall. This means that even if $p_zL = H$, we have $\dot{\mathbf{p}}^J \cdot \mathbf{n} = 0$ and then the
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Comparing the orientation development for two different degrees of confinement: H/L = 0.2 (top) and H/L = 0.3 (bottom): lubrication fully-confined macroscopic $\tilde{a}^C(t; H)$ versus confined Jeffery's $\mathbf{a}^{\text{disc}}(t)$.</figcaption></figure>

contact force vanishes. The rod having thus reached a permitted, unconfined Jeffery trajectory, it follows it for the remainder of the simulation. Although the final orientation, i.e. full alignment of the rod with the flow, is the same as the one reached by the purely unconfined Jeffery motion, the confined trajectory has a slightly higher time of flight.

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

### 3.1.2. Fully-confined orientation tensor

Here, the initial confined orientation is given by $\mathbf{p}^T(t = 0) = (-p_x, p_y, H/L)$, with $p_y \approx 0$, $H/L = 0.8$ and such that $\|\mathbf{p}(t = 0)\| = 1$. The confined Jeffery equation (14) was integrated with the initial condition $\mathbf{p}(t = 0)$. The orientation tensor at each time step

<table><caption>**Table 1**
Orientation time delay for different degrees of confinement $H/L$.</caption><thead><tr><td>$H/L$</td><td>$\Delta t / \Delta t_{fl}^j$</td></tr></thead><tbody><tr><td>0.8</td><td>1.3</td></tr><tr><td>0.6</td><td>1.9</td></tr><tr><td>0.4</td><td>2.6</td></tr><tr><td>0.2</td><td>4.1</td></tr></tbody></table>

$\mathbf{a}^{\text{disc}}(t)$ was calculated from $\mathbf{a}^{\text{disc}}(t) = \mathbf{p}(t) \otimes \mathbf{p}(t)$ (having all rods oriented in the same direction is equivalent to consider $N = 1$ in Eq. (17)), and compared with the fully-confined orientation tensor $\mathbf{a}^C(t)$ obtained by integrating Eq. (33) from the initial condition $\mathbf{a}^C(t = 0) = \mathbf{p}(t = 0) \otimes \mathbf{p}(t = 0)$ that allowed to define $\mathbf{b}(t = 0)$ and $\mathbf{c}(t = 0)$.

Fig. 4 compares both solutions. It can be noticed that they are in perfect agreement at the beginning of the orientation process when confinement is intense, while they differ slightly from each other at the end when fibres detach. As just mentioned, the confined macroscopic model involving $\mathbf{a}^C$ prevents detachment.

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

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

## 3.2. Evolution of a population of rods

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

### 3.2.1. Discrete calculation of the orientation tensor

We consider a population of $N = 4000$ rods having an almost isotropic initial orientation on the allowed part of the unit ball surface, the one limited by the parallels $z = \pm H$, as shown in Fig. 6. The trajectory of each $\mathbf{p}_i(t)$ is obtained by integrating the confined Jeffery model (14). The results will be compared with those of the unconfined case, $\mathbf{p}_i^j(t)$, governed by the standard Jeffery equation (1).

At each instant, the second-order orientation tensor is computed for both the confined and the unconfined systems, $\mathbf{a}^{\text{disc}}(t)$ and $\mathbf{a}^{j, \text{disc}}(t)$ respectively, according to:

$$ \mathbf{a}^{\text{disc}}(t) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{p}_i(t) \otimes \mathbf{p}_i(t), \qquad (42) $$

and

$$ \mathbf{a}^{j, \text{disc}}(t) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{p}_i^j(t) \otimes \mathbf{p}_i^j(t). \qquad (43) $$

In Fig. 7, one notices that the orientation $a_{xx}^{\text{disc}}$ experiences a slight delay with respect to $a_{xx}^{j, \text{disc}}$, and the peak in $a_{zz}^{j, \text{disc}}$ disappears when confinement effects act.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Fully-confined solutions $a^C(t)$ for different degrees of confinement $H/L$.</figcaption></figure>

When considering a population of rods rather than a single rod, the delaying effect of confinement is less noticeable due to the averaging of results involving a number of unconfined trajectories associated with rods that never reach the walls. Thus, we can conclude that the impact of confinement is quite moderate in what concerns the orientation moments.

### 3.2.2. Impact of the closure relation

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

The solution of the fully-confined macroscopic model for $a^C(t)$ is compared with the integration of the confined Jeffery model (14). Fig. 9 compares both solutions and differences are noticed from the very beginning of the orientation process. These differences could be attributed to the closure relations (38) involved in the formulation of $a^C$. It is easy to verify that the closure relations previously introduced for expressing $\mathbb{B}$ and $\mathbf{B}$ as a function of $b$ and $c$ are only exact in the case of full alignment (when the orientation pdf reduces to a Dirac distribution).

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

### 3.2.3. Unconfined versus confined Jeffery models

The results just discussed indicate that confinement only has a slight influence on the kinematics of the orientation process, as Fig. 7 reveals. It is thus of interest to have a further look at the un-confined Jeffery equation (1) and the associated evolution equation (20) for the orientation tensor $a^J$ when considering a confined initial orientation distribution. Here, we use the quadratic and hybrid closure relations for expressing the fourth-order orientation tensor $A^J$ in terms of $a^J$ [27]. Other closures exist, such as the natural and the orthotropic ones [10,11,13]. The former does not, however, have an explicit expression in the 3D case and the latter predicts spurious oscillations in absence of diffusion. For these reasons, we consider the quadratic are hybrid closures that are expected to perform reasonably well in unconfined conditions.

In the numerical experiments that follow, we specify initially two Gaussian orientation states on the allowed part of the unit ball surface, both depicted in Fig. 11, limited by the parallels $z = \pm H$,

with $H/L = 0.2$ and $0.8$. Thus, for integrating Eq. (20), we specify

$$ \mathbf{a}^J(t=0) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{p}_i(0) \otimes \mathbf{p}_i(0), \qquad (44) $$

where the initially-confined discrete orientations $\mathbf{p}_i(0)$ are those shown in Fig. 11.

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

* The out-of-plane component ($\bullet_{zz}$) for the unconfined Jeffery model is only slightly different from the one related to the confined model;

* For confined initial conditions, confined kinematics delay slightly the evolution of the second-order orientation tensor obtained from a population of rods large enough to consider the computed solution as almost exact statistically (i.e. solving the corresponding Fokker-Planck equation would give the same result). The delay increases with the degree of confinement, as measured by the reciprocal of the ratio $H/L$;

* Closure relations are responsible for artificially accelerating the evolution of the second-order moment resulting from the integration of its unconfined evolution equation with confined initial conditions. Closure approximations become less accurate (i.e. the orientation process accelerates) as the degree of confinement increases (i.e. as $H/L$ decreases);

* The quadratic and hybrid closures produce similar results.

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

A natural route is the development of empirical closure relations, inspired from the works of [10,11,13] that were successfully implemented in [23,29]. This could be a valuable route indeed, and a fitted closure should of course work well as long as the operating conditions remain similar to the ones that served to construct it. Although useful from an applicative point of view, this approach has two main limitations: (i) it could hide the real physical reasons for the noticed deviations; and (ii) a closure fitted from a particular confined flow could fail as soon as the flow and degree of confinement differ significantly from the ones that were used to obtain the particular fit.

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

Following the rationale described in [3], we attempted to describe the initial confined distribution depicted in Fig. 6 from its associated moments $\mathbf{a}$, $\mathbf{A}$ and $\mathbf{A}$. Use of the first three non-vanishing moments was found insufficient for approximating the
actual distribution. Many additional higher-order moments would be needed, which is not practical and would compromise the efficiency of the solution procedure. As discussed in [3], the orientation tensors can be viewed as the Fourier series expansion coefficients of the orientation distribution function. Limitations of Fourier series for approximating rectangular functions like the ones associated to a uniform distribution in a region of the conformation space (surface of unit sphere) are well known. In that situation, it is preferable to consider the direct solution of the Fokker-Planck equation, which constitutes an appealing route for further developments [8,9].

### 3.3. Lubrication approach

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

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

## 4. Conclusions

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

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

In view of this result, we performed the same simulations with the unconfined Jeffery equation (1) and the associated macroscopic evolution equation (20) for the second-order orientation tensor $\mathbf{a}^J$.

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

Thus, we conclude from this study that the main challenge with traditional macroscopic models involving moments of the orientation pdf lies more with representation capabilities in highly confined conditions than with a suitable description of the induced orientation kinematics. Use of the averaged fully-confined macroscopic descriptor $\tilde{\mathbf{a}}^C$ proposed in this paper is recommended in future theoretical developments.

## Acknowledgements

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

## Appendix A. Advanced mesoscopic modelling

For representing accurately the orientation distribution at the mesoscopic scale, we could consider two pdf's, $\psi^J(\mathbf{x}, t, \mathbf{p} \in \mathcal{J})$ and $\psi^C(\mathbf{x}, t, \mathbf{p} \in \mathcal{C})$, with

$$ \begin{cases} \mathcal{J} = \{\mathbf{p}, p_z \in (-H/L, H/L)\} \\ \mathcal{C} = \{\mathbf{p}, p_z = \pm H/L\} \end{cases}, \qquad (A.1) $$

where the normalization condition reads

$$ \int_{\mathcal{J}} \psi^J(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p} + \int_{\mathcal{C}} \psi^C(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p} = 1. \qquad (A.2) $$

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

$$ \left\{ \begin{aligned} \frac{\partial \psi^J}{\partial t} + \nabla_p \cdot (\dot{\mathbf{p}}^J \psi^J) &= 0, & \mathbf{p} \in \mathcal{J} \\ \frac{\partial \psi^C}{\partial t} + \nabla_p \cdot (\dot{\mathbf{p}}^M \psi^C) &= \mathcal{Q}^+ - \mathcal{Q}^-, & \mathbf{p} \in \mathcal{C}^+ \end{aligned} \right. \qquad (A.3) $$

where, due to symmetry considerations, we only consider the upper parallel $\mathcal{C}^+ = \{\mathbf{p}, p_z = H/L\}$, and $\dot{\mathbf{p}}^M$ represents the velocity on the manifold, defined from

$$ \dot{\mathbf{p}}^M = \begin{cases} \dot{\mathbf{p}} & \text{if } \dot{\mathbf{p}} \cdot \mathbf{n} = 0 \\ \mathbf{0} & \text{if } \dot{\mathbf{p}} \cdot \mathbf{n} \neq 0 \end{cases}, \qquad (A.4) $$

with $\dot{\mathbf{p}}$ given by Eq. (14).

Here, $\mathcal{Q}^+ = (\dot{\mathbf{p}}^+ \cdot \mathbf{t}) \psi^J(\mathbf{p})$ represents the unconfined rods reaching the manifold $\mathcal{C}^+$, $\mathbf{t}$ being the unit vector tangent to the unit sphere $\mathcal{S}$, normal to the manifold $\mathcal{C}^+$ and pointing outward of the allowed region $\mathcal{J}$, with the upstream velocity $\dot{\mathbf{p}}^+$ given by

$$ \dot{\mathbf{p}}^+ = \begin{cases} \dot{\mathbf{p}}^J & \text{if } \dot{\mathbf{p}}^J \cdot \mathbf{n} > 0 \\ \mathbf{0} & \text{otherwise} \end{cases}. \qquad (A.5) $$

On the other hand, confined rods leaving the manifold $\mathcal{C}^+$ are given by $\mathcal{Q}^-$, with $\mathcal{Q}^- = -(\dot{\mathbf{p}}^- \cdot \mathbf{t}) \psi^C(\mathbf{p})$, with

$$ \dot{\mathbf{p}}^- = \begin{cases} \dot{\mathbf{p}}^J & \text{if } \dot{\mathbf{p}}^J \cdot \mathbf{n} < 0 \\ \mathbf{0} & \text{otherwise} \end{cases}. \qquad (A.6) $$

The only boundary condition to be prescribed at the boundary of $\mathcal{J}$ to ensure conservation of probability reads

$$ \psi^J(\mathbf{p})|_{\partial \mathcal{J}^-} = \psi^C(\mathbf{p})|_{\partial \mathcal{J}^-}, \qquad (A.7) $$
where $\partial \mathcal{J}^-$ denotes the part of the boundary of $\mathcal{J}$ through which rods leaving the manifold $\mathcal{C}$ come into domain $\mathcal{J}$.

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

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

## Appendix B. Evolution equation for the confined orientation tensor

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

$$ \dot{\mathbf{p}} = \dot{\mathbf{p}}^J - \frac{1}{1-p_z^2} [\dot{\mathbf{p}}^J]_z (\mathbf{n} - p_z \mathbf{p}), \quad (B.1) $$

where $\mathbf{p}$ is written as

$$ \mathbf{p} = \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix}, \quad (B.2) $$

and the gradient of velocity as

$$ \nabla \mathbf{v} = \begin{pmatrix} \mathbf{G} & \mathbf{g} \\ \mathbf{j}^T & \mathcal{G} \end{pmatrix}. \quad (B.3) $$

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

$$ \nabla \mathbf{v} \cdot \mathbf{p} = \begin{pmatrix} \mathbf{G} & \mathbf{g} \\ \mathbf{j}^T & \mathcal{G} \end{pmatrix} \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix} = \begin{pmatrix} \mathbf{G} \cdot \mathbf{q} + \mathcal{G} p_z \\ \mathbf{j}^T \cdot \mathbf{q} + \mathcal{G} p_z \end{pmatrix}, \quad (B.4) $$

and

$$ \begin{aligned} (\mathbf{p}^T \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \mathbf{p} &= \begin{pmatrix} (\mathbf{q}^T & p_z) \begin{pmatrix} \mathbf{G} & \mathbf{g} \\ \mathbf{j}^T & \mathcal{G} \end{pmatrix} \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix} \end{pmatrix} \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix} \\ &= (\mathbf{q}^T \cdot \mathbf{G} \cdot \mathbf{q} + (\mathbf{q}^T \cdot \mathbf{g}) p_z + (\mathbf{j}^T \cdot \mathbf{q}) p_z + \mathcal{G} p_z^2) \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix}. \end{aligned} \quad (B.5) $$

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

$$ [\dot{\mathbf{p}}^J]_z = \mathbf{j}^T \cdot \mathbf{q} + \mathcal{G} p_z - (\mathbf{q}^T \cdot \mathbf{G} \cdot \mathbf{q}) p_z - (\mathbf{q}^T \cdot \mathbf{g}) p_z^2 - (\mathbf{j}^T \cdot \mathbf{q}) p_z^2 - \mathcal{G} p_z^3 \quad (B.6) $$

and

$$ (\mathbf{n} - p_z \mathbf{p}) = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix} - p_z \begin{pmatrix} \mathbf{q} \\ p_z \end{pmatrix} = \begin{pmatrix} -p_z \mathbf{q} \\ 1 - p_z^2 \end{pmatrix}. \quad (B.7) $$

Thus, we finally obtain

$$ \begin{aligned} \dot{\mathbf{q}} = \mathbf{G} \cdot \mathbf{q} + \mathcal{G} p_z - (\mathbf{q}^T \cdot \mathbf{G} \cdot \mathbf{q}) \mathbf{q} - (\mathbf{q}^T \cdot \mathbf{g}) p_z \mathbf{q} - (\mathbf{j}^T \cdot \mathbf{q}) p_z \mathbf{q} - \mathcal{G} p_z^2 \mathbf{q} \\ - \frac{1}{1-p_z^2} \left( -(\mathbf{j}^T \cdot \mathbf{q}) p_z \mathbf{q} - \mathcal{G} p_z^2 \mathbf{q} + (\mathbf{q}^T \cdot \mathbf{G} \cdot \mathbf{q}) p_z^2 \mathbf{q} \right. \\ \left. + (\mathbf{q}^T \cdot \mathbf{g}) p_z^3 \mathbf{q} + (\mathbf{j}^T \cdot \mathbf{q}) p_z^3 \mathbf{q} + \mathcal{G} p_z^4 \mathbf{q} \right), \end{aligned} \quad (B.8) $$

that can be rewritten as

$$ \dot{\mathbf{q}} = \mathbf{G} \cdot \mathbf{q} - \delta_1 (\mathbf{q}^T \cdot \mathbf{G} \cdot \mathbf{q}) \mathbf{q} + \delta_2 \mathbf{g} - \delta_3 (\mathbf{q}^T \cdot \mathbf{g}) \mathbf{q} \quad (B.9) $$

with $\delta_1 = \frac{1}{1-p_z^2}$, $\delta_2 = p_z$ and $\delta_3 = \frac{p_z}{1-p_z^2}$.

## Appendix C. Time evolution of the averaged confined orientation tensor based on the lubrication approximation

By introducing Eqs. (37) and (40) into the time derivative of Eq. (41), we obtain

$$ \begin{aligned} \dot{\tilde{\mathbf{b}}} = \mathbf{G} \cdot \tilde{\mathbf{b}} + \tilde{\mathbf{b}} \cdot \mathbf{G}^T - 2\mathbf{G} : \int_0^H \omega(z) \delta_1 \mathbf{B} \, dz \\ + (\mathbf{g} \otimes \tilde{\mathbf{c}} + \tilde{\mathbf{c}} \otimes \mathbf{g}) - 2 \left( \int_0^H \omega(z) \delta_3 \mathbf{B} \, dz \right) \cdot \mathbf{g}, \end{aligned} \quad (C.1) $$

and

$$ \begin{aligned} \dot{\tilde{\mathbf{c}}} = \mathbf{G} \cdot \tilde{\mathbf{c}} - \mathbf{G} : \int_0^H \omega(z) \delta_1 \delta_2 \mathbf{B} \, dz \\ + \left( \int_0^H \omega(z) \delta_2^2 \, dz \right) \mathbf{g} - \left( \int_0^H \omega(z) \delta_3 \delta_2 \mathbf{b} \, dz \right) \cdot \mathbf{g}, \end{aligned} \quad (C.2) $$

with $p_z = z/L$ and $\delta_1 = \frac{1}{1-p_z^2}$, $\delta_2 = p_z$ and $\delta_3 = \frac{p_z}{1-p_z^2} = \delta_1 \delta_2$.

A closed solution cannot be derived, but approximate expressions can be obtained by closing the integral terms. For that purpose, we define

$$ \left\{ \begin{aligned} \tilde{\mathbf{b}} &= \int_0^H \omega(z) \mathbf{b} \, dz \\ \tilde{\mathbf{c}} &= \int_0^H \omega(z) \mathbf{c} \, dz \\ \tilde{\mathbf{B}} &= \int_0^H \omega(z) \mathbf{B} \, dz \\ \tilde{\mathbf{B}} &= \int_0^H \omega(z) \mathbf{B} \, dz \end{aligned} \right. , \quad (C.3) $$

leading to the decomposition

$$ \left\{ \begin{aligned} \mathbf{b} &= \tilde{\mathbf{b}} + \Delta\mathbf{b} \\ \mathbf{c} &= \tilde{\mathbf{c}} + \Delta\mathbf{c} \\ \mathbf{B} &= \tilde{\mathbf{B}} + \Delta\mathbf{B} \\ \mathbf{B} &= \tilde{\mathbf{B}} + \Delta\mathbf{B} \end{aligned} \right. , \quad (C.4) $$

with, by construction,

$$ \left\{ \begin{aligned} \int_0^H \omega(z) \Delta\mathbf{b} \, dz &= 0 \\ \int_0^H \omega(z) \Delta\mathbf{c} \, dz &= 0 \\ \int_0^H \omega(z) \Delta\mathbf{B} \, dz &= 0 \\ \int_0^H \omega(z) \Delta\mathbf{B} \, dz &= 0 \end{aligned} \right. , \quad (C.5) $$

The co-factors involving the delta coefficients in Eqs. (C.1) and (C.2) are noted for the sake of notational simplicity as $\mathcal{D}_1 = \delta_1$, $\mathcal{D}_2 = \delta_3$, $\mathcal{D}_3 = \delta_1 \delta_2$, $\mathcal{D}_4 = \delta_2^2$ and $\mathcal{D}_5 = \delta_2 \delta_3$. They accept the decomposition

$$ \mathcal{D}_i = \tilde{\mathcal{D}}_i + \Delta\mathcal{D}_i, \quad i = 1, \dots, 5, \quad (C.6) $$

with $\int_0^H \omega(z) \Delta\mathcal{D}_i \, dz = 0, \forall i$.

Thus, neglecting integrals involving products of variations, as for example $\int_0^H \omega(z) \Delta\mathcal{D}_1 \Delta b \, dz \approx 0$, and taking into account that tilda variables do not depend on the z-coordinate and that $\int_0^H \omega(z) \, dz = 1$, Eqs. (C.1) and (C.2) read:
$$ \dot{\tilde{\mathbf{b}}} = \mathbf{G} \cdot \tilde{\mathbf{b}} + \tilde{\mathbf{b}} \cdot \mathbf{G}^T - 2\tilde{\mathbf{D}}_1 \mathbf{G} : \tilde{\mathbf{B}} + (\mathbf{g} \otimes \tilde{\mathbf{c}} + \tilde{\mathbf{c}} \otimes \mathbf{g}) - 2\tilde{\mathbf{D}}_2 \tilde{\mathbf{B}} \cdot \mathbf{g}. \quad (C.7) $$

and

$$ \dot{\tilde{\mathbf{c}}} = \mathbf{G} \cdot \tilde{\mathbf{c}} - \tilde{\mathbf{D}}_3 \mathbf{G} : \tilde{\mathbf{B}} + \tilde{\mathbf{D}}_4 \mathbf{g} - \tilde{\mathbf{D}}_5 \tilde{\mathbf{b}} \cdot \mathbf{g}. \quad (C.8) $$

These equations require appropriate closure relations for the average higher-order orientation moments $\tilde{\mathbf{B}}$ and $\tilde{\mathbf{B}}$.

## References

[1] E. Abisset-Chavanne, F. Chinesta, J. Ferec, G. Ausias, R. Keunings, On the multiscale description of dilute suspensions of non-brownian rigid clusters composed of rods, Journal of Non-Newtonian Fluid Mech. 222 (2015) 34–44.

[2] E. Abisset-Chavanne, J. Ferec, G. Ausias, E. Cueto, F. Chinesta, R. Keunings, A second-gradient theory of dilute suspensions of flexible rods in a newtonian fluid, Arch. Comput. Methods Eng. 22 (2015) 511–527.

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

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

[5] C. Binetruy, F. Chinesta, R. Keunings, Flows in Polymers, Reinforced Polymers and Composites: A Multiscale Approach, Springer, Springerbriefs, 2015.

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

[7] D. Borzacchiello, E. Abisset-Chavanne, F. Chinesta, R. Keunings, Orientation kinematics of short fibres in a second-order viscoelastic fluid, Rheol. Acta (submitted for publication),

[8] F. Chinesta, A. Ammar, A. Leygue, R. Keunings, An overview of the proper generalized decomposition with applications in computational rheology, J. Non-Newtonian Fluid Mech. 166 (2011) 578–592.

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

[10] D.H. Chung, T.H. Won, Improved model of orthotropic closure approximation for flow induced fiber orientation, Polym. Compos. 22/5 (2001) 636–649.

[11] J.S. Cintra, Ch. L. Tucker III, Orthotropic closure approximations for flow-induced fiber orientation, J. Rheol. 39/6 (1995) 1095–1121.

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

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

[14] J. Ferec, G. Ausias, M.C. Heuzey, P. Carreau, Modeling fiber interactions in semi-concentrated fiber suspensions, J. Rheol. 53/1 (2009) 49–72.

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

[16] E. Gavze, E. Shapiro, Particles in a shear flow near a solid wall: effects of non sphericity on forces and velocities, Int. J. Multiph. Flow 23/1 (1997) 155–182.

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

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

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

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

[21] C. Jayageeth, V.I. Sharma, A. Singh, Dynamics of short fiber suspensions in bounded shear flow, Int. J. Multiph. Flow 35 (2009) 261–269.

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

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

[24] T.-H. Le, P. Dumont, L. Orgeas, D. Favier, L. Salvo, E. Boller, X-ray phase contrast microtomography for the analysis of the fibrous microstructure of SMC composites, Compos.: Part A 39 (2008) 91–103.

[25] M.B. Mackaplowt, E.S.G. Shaqfeh, A numerical study of the rheological properties of suspensions of rigid, non-brownian fibres, J. Fluid Mech. 329 (1996) 155–186.

[26] A. Ozolins, U. Strautins, Simple models for wall effect in fiber suspension flows, Math. Model. Anal. 19/1 (2014) 75–84.

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

[28] J. Phelps, Ch. Tucker, An anisotropic rotary diffusion model for fiber orientation in short and long fiber thermoplastics, J. Non-Newtonian Fluid Mech. 156/3 (2009) 165–176.

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

[30] C.A. Stover, C. Cohen, The motion of rodlike particles in the pressure-driven flow between two flat plates, Rheol. Acta 29 (1990) 192–203.

[31] J. Wang, C.A. Silva, J.C. Viana, F.W.J. van Hattum, A.M. Cunha, Ch. Tucker, Prediction of fiber orientation in a rotating compressing and expanding mold, Polym. Eng. Sci. (2008) 1405–1413.

[32] J. Wang, J. O’Gara, Ch. Tucker, An objective model for slow orientation kinetics in concentrated fiber suspensions: theory and rheological evidence, J. Rheol. 52/5 (2008) 1179–1200.