# A Second-Gradient Theory of Dilute Suspensions of Flexible Rods in a Newtonian Fluid

Emmanuelle Abisset-Chavanne, Julien Férec, Gilles Ausias, Elías Cueto, Francisco Chinesta, Roland Keunings

► To cite this version:

Emmanuelle Abisset-Chavanne, Julien Férec, Gilles Ausias, Elías Cueto, Francisco Chinesta, et al. A Second-Gradient Theory of Dilute Suspensions of Flexible Rods in a Newtonian Fluid. Archives of Computational Methods in Engineering, 2015, 22 (3), pp.511-527. 10.1007/s11831-014-9128-6. hal-01515081

HAL Id: hal-01515081

https://hal.science/hal-01515081v1

Submitted on 27 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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Second-Gradient Theory of Dilute Suspensions of Flexible Rods in a Newtonian Fluid

E. Abisset-Chavanne¹, J. Férec², G. Ausias²,
E. Cueto³, F. Chinesta¹*, R. Keunings⁴

**Abstract** Most suspension descriptions used nowadays are based on Jeffery's model or some phenomenological modifications of it that do not take into account size effects: the kinematics and stresses do not involve a micro-mechanical characteristic length and thus, the predicted rheological properties are necessarily independent of rod length. More sophisticated models able to enrich first-gradient kinematics as well as to activate rod-bending mechanisms are needed, in particular to explain the mild elasticity observed experimentally. In this paper we propose a second-gradient theory for dilute suspensions of flexible rods in a Newtonian fluid that is indeed able to activate rod bending and predict viscoelastic behaviour.

**Keywords:** Rod suspensions, Jeffery's equation, Size effects, Second gradient, Rod bending, hyper-stress, viscoelasticity

Emmanuelle Abisset-Chavanne & Francisco Chinesta
GeM, UMR CNRS-Centrale Nantes
1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France
E-mail: Emmanuelle.Abisset-Chavanne@ec-nantes.fr, Francisco.Chinesta@ec-nantes.fr

Julien Ferec & Gilles Ausias
LIMATB, Université de Bretagne-Sud
rue de St Maudé, F-56325 Lorient, France
E-mail: julien.ferec@univ-ubs.fr, gilles.ausias@univ-ubs.fr

Elias Cueto
I3A, Aragon Institute of Engineering Research
Universidad de Zaragoza, Maria de Luna s/n, Zaragoza, Spain
E-mail: Elias.Cueto@unizar.es

Roland Keunings
ICTEAM, Université catholique de Louvain
Bat. Euler, Av. Georges Lemaitre 4, B-1348 Louvain-la-Neuve, Belgium
E-mail: roland.keunings@uclouvain.be

<footnote>* Francisco Chinesta, corresponding author</footnote>
# 1 Introduction

Short fibers are widely used for reinforcing polymers. The resulting microstructure has been traditionally described at three different scales.

At the finest scale the microstructure can be described, when the suspension is dilute enough, by tracking a population of rods that move with the suspending fluid (assumed Newtonian) and orient depending on the velocity gradient according to Jeffery's equation [22] that relates the orientation evolution with the flow velocity field, in particular to its first gradient. In that case the motion and orientation of each fiber is assumed decoupled from the others. When the concentration increases rod-rod interactions cannot be neglected anylonger. They must be included into the force and moment balances of each rod in order to describe their motions.

When the representative population involves too many fibers, the computational efforts to track the population becomes unaffordable. Thus, the simple and well-defined physics must be sacrificed in order to derive coarser descriptions. Kinetic theory approaches [6] [11] [24] describe such systems at the mesoscopic scale. Their main advantage is their ability to address macroscopic systems, while keeping track of 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 physical space, time and a number of conformational coordinates, i.e. rod orientation in the case of slender bodies suspensions.

The moments of the orientation distribution function constitute a coarser description in general used in macroscopic modeling [1]. At the macroscopic scale the equations governing the time evolution of these moments usually involve closure approximations whose impact on the results is unpredictable [7] [12] [23]. Alternatively, macroscopic equations are carefully postulated in order to guarantee the model objectivity and its thermodynamical admissibility.

In the case of dilute suspensions of short fibers, the three scales have been extensively considered to model the associated systems without major difficulties. However, as soon as the concentration increases, difficulties appear. In the semi-dilute and semi-concentrated regimes fiber-fiber interactions occur, but in general they can be accurately modeled by introducing a sort of randomizing diffusion term [15]. There is a wide literature concerning dilute and semi-dilute suspensions, addressing modeling [5] [17] [18] [19] [20], flows [3] [4] [8] [35] and rheology [27] [29] [30].

Thus, most of suspension descriptions used nowadays are based on Jeffery's model or some phenomenological modification of it that do not take into account size effects: the kinematics and stresses do not involve a micro-mechanical characteristic length and thus, the predicted rheological properties are necessarily independent of rod length.

Size effects disappear as soon as a constant gradient of the fluid velocity is postulated at the scale of the rod (first-gradient framework). In that case the rod kinematics consist of an affine transformation corrected in order to ensure rod inextensibility. Obviously, because of the transformation linearity, rod bending mechanisms are prevented. In [21], the deformation of an inextensible thread without bending rigidity was considered, and it was shown that in a flow with constant velocity gradient, deformations are activated as soon as a nearly straight initial configuration is assumed. In order to activate bending mechanisms, we considered in our former works first-gradient hydrodynamics acting on rods with non-straight unloaded configurations [9] [10]. In [34], the
authors considered a three-bead-two-rod system equipped with an elastic spring that resists the rod differential rotation.

In this paper, we propose a second-gradient description of rod suspensions that allows activation of rod bending and thus naturally introduces elastic stresses, although it considers straight rods and predicts Jeffery rod kinematics.

In the next section, we develop first and second-gradient models for the kinematics of rods immersed in a Newtonian fluid, with the rods being modeled as dumbbells. We will prove that in both descriptions the rod kinematics are described by Jeffery's equation and that, with forces acting along the rod direction, rod bending cannot be activated. In section 3, we consider an alternative approach in which hydrodynamic forces are distributed all along the rod length. In that case, the rod kinematics are again described by Jeffery's equation but, as soon as the second gradient is considered, rod bending is activated. From the microscopic description, we derive in Section 4 a second-gradient macroscopic model that couples the rod behaviour with the flow kinematics. In section 5, we show that the proposed theory predicts a Maxwell linear viscoelastic behaviour. Section 6, we study the model predictions in the inception and cessation of simple shear flow. Finally, in the spirit of [25], we briefly address in Section 7 the complex flow in a driven cavity assuming coupling between rod and flow kinematics.

We voluntarily omit a detailed numerical analysis of the proposed theory in order to limit the size of this paper. It will be reported in future publications. We thus mainly focus our analysis on the effect of flow kinematics on the rod conformation, that constitutes the main originality of our proposal.

*Remark 1.* In this paper, we consider the following tensor products:

- if **a** and **b** are first-order tensors then the single contraction “·” reads (a·b) = a_j b_j (Einstein's summation convention);

- 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 first-order tensors then the cross product “×” reads (a × b)_j = ε_{jmn} a_m b_n (Einstein's summation convention) with ε_{jmn} the components of the Levi-Civita tensor **ε** (also known as permutation tensor);

- if **a** and **b** are respectively second and first-order tensors then the single contraction “·” reads (a·b)_j = a_j b_m (Einstein's summation convention);

- if **a** and **b** are second-order tensors then the single contraction “·” reads (a·b)_{jk} = a_j b_m (Einstein's summation convention);

- if **a** and **b** are second-order tensors then the double contraction “::” reads (a : b) = a_jk b_{kj} (Einstein's summation convention);

- if **a** and **b** are third-order tensors then the triple contraction “.“·” reads (a···b) = a_jkm b_{mkj} (Einstein's summation convention);

- if **a** and **b** are fourth-order tensors then the fourth contraction “’::” reads (a::b) = a_jkmn b_{nmkj} (Einstein's summation convention).

## 2 From first to second-gradient descriptions of rod kinematics

We consider a suspending medium consisting of a Newtonian fluid of viscosity η in which there are suspended rigid and non-Brownian rods. We assume as first approximation that their presence and orientation do not affect the flow kinematics (this hypoth-
<figure><img src="image_2.png" /><figcaption>Fig. 1 Hydrodynamic forces applied on a rod immersed in a Newtonian fluid</figcaption></figure>

esis will be relaxed later) that is defined by the velocity field $v(\mathbf{x}, t)$, with $\mathbf{x} \in \Omega \in \mathbb{R}^d$, $d = 2, 3$.

The conformation of each rod of length $2L$ can be described from its orientation, given by the unit vector $\mathbf{p}$ located at the rod center of gravity $G$ and aligned along the rod axis. The rod's mass is neglected and with its micro-scale inertial effects.

## 2.1 First-gradient modeling: Jeffery's equation

The orientation evolution of an ellipsoidal particle suspended in a Newtonian fluid is described by Jeffery's equation [22], that for rods (ellipsoids with infinite aspect ratio) can be derived by considering the system illustrated in Fig. 1 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}$ depend on the difference of velocities between the fluid and the bead, the first one given by $\mathbf{v}_0 + \nabla \mathbf{v} \cdot \mathbf{p}L$ and the second one by $\mathbf{v}_G + \dot{\mathbf{p}}L$. Thus, force $\mathbf{F}(\mathbf{p}L)$ reads:

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

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

Obviously if $\mathbf{F}$ applies on the bead $\mathbf{p}L$, then at the opposite bead $-\mathbf{p}L$ the resulting force reads

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

By adding Eqs. (1) and (2) and enforcing the force balance (neglecting inertia effects, as the rod mass is assumed negligible), we obtain

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

than implies $\mathbf{v}_0 = \mathbf{v}_G$, that is, the rod center of gravity is moving with the fluid velocity. For alleviating the notation, we consider $\mathbf{F} = \mathbf{F}(\mathbf{p}L)$ and $\mathbf{F}(-\mathbf{p}L) = -\mathbf{F}$.
As the resulting torque must also vanish, the only possibility is that force **F** acts along the direction defined by **p**, that is **F** = $\lambda\mathbf{p}$, with $\lambda \in \mathbb{R}$. Thus we can write

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

Premultiplying Eq. (4) by **p** and taking into account that **p**·**p** = 1 and consequently **p**·**p** = 0, we get:

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

that allows us to write

$$ \xi L (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p} = \xi L (\nabla \mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}), \quad (6) $$

from which we finally obtain the classical Jeffery equation

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

*Remark 2.* Because the factor $\xi L$ appears in both sides of Eq. (6), the predicted rod kinematics does not involve size effects.

The forces applied at the rod ends **p**L and $-\mathbf{p}L$ are respectively $\lambda\mathbf{p}$ and $-\lambda\mathbf{p}$, both directed along the rod direction and by construction auto-equilibrated.

With $\lambda$ given by Eq. (5) the force **F** reads

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

which is aligned, as expected, in the rod direction.

Such a force generates the rod tension or compression depending on the sign of $\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})$. In the case of a flexible rod it can produce its extension or compression, the last being able to produce buckling.

Thus, a certain elasticity could be expected, however for standard materials the elasticity related to the axial tension is too small and compression mechanics leads to buckling and the consequent rod rupture.

Rod bending seems to be another plausible mechanism responsible for elastic effects, however forces aligned in the rod direction cannot activate rod bending. The introduction of such effects requires at least a second-gradient theory, as described in what follows.

## 2.2 Second-gradient description with concentrated forces at the rod beads

We now consider the same system but with a higher-order description of the fluid velocity field at the rod scale. Again, forces applied on each bead **F** depend on the difference of velocities between the fluid and the bead, the first one now including the second-order velocity gradient **H** according to $\mathbf{v}(\mathbf{p}L) = \mathbf{v}_0 + \nabla \mathbf{v} \cdot \mathbf{p}L + (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p})) L^2$ and the second one given by $\mathbf{v}_G + \dot{\mathbf{p}}L$. Thus, the force $\mathbf{F}(\mathbf{p}L)$ reads:

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

where the third-order tensor **H** is defined by $\mathbf{H}_{ijk} = \frac{1}{2} \frac{\partial^2 \mathbf{v}_i}{\partial x_j \partial x_k}$.
Obviously, if **F** applies on the bead **pL**, then at the opposite bead **-pL** the resulting force reads

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

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

$$ \mathbf{F}(\mathbf{p}L) + \mathbf{F}(-\mathbf{p}L) = \mathbf{0}, \quad (11) $$

than implies $\mathbf{v}_0 - \mathbf{v}_G = -(\mathbf{H} : (\mathbf{p} \otimes \mathbf{p}))L^2$, that is, the rod center of gravity has a relative velocity with respect to the fluid at this position.

*Remark 3.* The fact of having obtained a non-zero relative velocity $\mathbf{v}_0 - \mathbf{v}_G \neq \mathbf{0}$ does not imply the existence of a migration mechanism, as shown in Appendix A. See also [33].

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

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

which yields the same Jeffery equation that was obtained when considering the first-order velocity gradient:

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

The forces being again aligned in the rod direction, one could infer that the second gradient does not suffice for activating bending mechanisms.

Consideration of a third-gradient description predicts rod kinematics that differ from those dictated by Jeffery's equation (see Appendix B). However, the forces remain aligned along the rod direction and rod bending thus cannot be activated either.

Until now, forces were assumed applied at the rod ends (beads). However, we can distribute them all along the rod length as commonly considered when calculating particles motion by using DPD (dissipative particle dynamics) methods.

In the next section, we consider forces applied all along the rod length as was considered in [32]. We prove that as soon as a second-gradient description is retained, bending mechanism appears naturally.

## 3 Second-gradient description of rods with distributed forces

We consider now the system illustrated in Fig. 2 consisting of a rod and the hydrodynamical forces applied all along its length. With the same reasoning and notation as above, the applied distributed force $\mathbf{f}(s)$ at position $s\mathbf{p}$, $s \in [-L, L]$ reads

$$ \mathbf{f}(s) = \bar{\xi} (\mathbf{v}_0 + \nabla\mathbf{v} \cdot \mathbf{p}s + (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p}))s^2 - \mathbf{v}_G - \dot{\mathbf{p}}s). \quad (14) $$

The resultant force **F** must vanish, that is

$$ \mathbf{F} = \int_{-L}^{+L} \mathbf{f}(s) \, ds = \mathbf{0}, \quad (15) $$
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Distributed hydrodynamic forces applied on a rod immersed in a Newtonian fluid, considering the second-order velocity gradient.</figcaption></figure>

that leads to the following expression for the sliding velocity (rod-fluid relative velocity at the rod center of gravity):

$$ \mathbf{v}_0 - \mathbf{v}_G = -\frac{L^2}{3} (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p})). \quad (16) $$

Thus, the distributed force reads:

$$ \mathbf{f}(s) = \bar{\xi} \left( \nabla \mathbf{v} \cdot \mathbf{p}s + (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p})) \left( s^2 - \frac{L^2}{3} \right) - \dot{\mathbf{p}}s \right), \quad (17) $$

which leads to the moment $\mathbf{m}(s)$

$$ \mathbf{m}(s) = s\mathbf{p} \times \bar{\xi} \left( \nabla \mathbf{v} \cdot \mathbf{p}s + (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p})) \left( s^2 - \frac{L^2}{3} \right) - \dot{\mathbf{p}}s \right), \quad (18) $$

from which we can evaluate the resultant moment $\mathbf{M}$

$$ \mathbf{M} = \int_{-L}^{+L} \mathbf{m}(s) \, ds = \mathbf{p} \times \bar{\xi} \, \alpha (\nabla \mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}) = \mathbf{0}, \quad (19) $$

with $\alpha = \int_{-L}^{+L} s^2 \, ds = \frac{2}{3} L^3$. This again yields Jeffery's equation

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

i.e. the same equation as that obtained by assuming forces applied at the rod beads.

Introducing Jeffery's equation (20) into the distributed force expression (17), we obtain

$$ \mathbf{f}(s) = \bar{\xi} \left( (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p})) \left( s^2 - \frac{L^2}{3} \right) + (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}s \right). \quad (21) $$

This force can be decomposed into two components, i.e. one, $\mathbf{f}^{\parallel}(s)$, aligned with the rod and the other, $\mathbf{f}^{\perp}(s)$, perpendicular to it:

$$ \mathbf{f}^{\parallel}(s) = \bar{\xi} \left( (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p})) \mathbf{p} \left( s^2 - \frac{L^2}{3} \right) + (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}s \right), \quad (22) $$

and

$$ \mathbf{f}^{\perp}(s) = \bar{\xi} \left[ (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p})) - (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p})) \mathbf{p} \right] \left( s^2 - \frac{L^2}{3} \right). \quad (23) $$
Remark 4. We can notice from Eq. (23) that, when considering distributed forces within a first-gradient description ($\mathbf{H} = \mathbf{0}$), the perpendicular component of the resulting distributed forces vanishes, i.e. $\mathbf{f}^{\perp}(s) = \mathbf{0}$ and bending mechanisms are once again absent. However, bending seems possible as soon as second-gradient descriptions implying $\mathbf{H} \neq \mathbf{0}$ are retained.

The resultant axial force $\mathbf{F}^{\parallel}$ reads

$$ \mathbf{F}^{\parallel} = \int_{0}^{L} \mathbf{f}^{\parallel}(s) \, ds = \bar{\xi} \frac{L^2}{2} (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p}, \quad (24) $$

expression that corresponds to the one obtained in the case of concentrated forces if we consider the following relation between the distributed and concentrated friction coefficients:

$$ \xi = \frac{L}{2} \bar{\xi}. \quad (25) $$

The main difference with respect to the situation in which forces were assumed applied at the rod beads, occurs when considering the distributed force perpendicular to the rod. In this case

$$ \mathbf{F}^{\perp} = \int_{0}^{L} \mathbf{f}^{\perp}(s) \, ds = \mathbf{0}, \quad (26) $$

in agreement with the results found when considering concentrated forces, but the distributed force implies the existence of a bending moment $\mathcal{M}(s)\mathbf{k}$ (acting in the out-of-plane direction defined by the unit vector $\mathbf{k}$)

$$ \mathcal{M}(s)\mathbf{k} = \int_{s}^{L} (r-s)\mathbf{p} \times \mathbf{f}^{\perp}(r) \, dr. \quad (27) $$

When the rod is rigid there are no major consequences, but in the case of flexible rods this force creates rod bending with the associated elastic effects.

Remark 5. In the case of elastic rods, the bending moment implies rod curvature. Within the small strain and displacement hypotheses, the curvature is given by

$$ \frac{d^2 u}{ds^2} = \frac{\mathcal{M}(s)}{EI}, \quad (28) $$

where $u(s)$ is the rod deflection with respect to its undeformed configuration, $E$ the elastic modulus and $I$ the moment of area of the rod cross section with respect to the out-of-plane direction. By integrating this equation twice, we can obtain the rod bent configuration that, as expected, is symmetric with respect to the rod center of gravity.

It is important to notice that moments imply an extra length $L$ with respect to forces, and moreover each integration step for moving from curvature to displacement introduces the length $L$ as factor.

Having proved that second-gradient descriptions activate bending effects, we propose in the next section a theoretical model able to couple flow and rod kinematics within a second-gradient framework.
# 4 Second-gradient flow model of dilute suspensions involving flexible rods

## 4.1 Second-gradient fluid description

Within the first-gradient formulation, the internal power for a Newtonian fluid $W_{int}$ reads

$$W_{int} = \int_{\Omega} \mathbf{T} : \nabla \mathbf{v} \, d\mathbf{x}, \qquad (29)$$

where $\mathbf{T} = -p\mathbf{I} + \boldsymbol{\tau}$ with $\boldsymbol{\tau} = 2\eta\mathbf{D}$. The associated balance of momentum reads:

$$\rho \dot{\mathbf{v}} = \nabla \cdot \mathbf{T}. \qquad (30)$$

Fried and Gurtin [16] proposed a second-gradient formulation involving the vorticity gradient. They proposed a non-standard form of the principle of virtual power with the internal power given by

$$W_{int} = \int_{\Omega} \mathbf{T} : \nabla \mathbf{v} \, d\mathbf{x} + \int_{\Omega} \mathbf{G} : \nabla \boldsymbol{\omega} \, d\mathbf{x}, \qquad (31)$$

where $\boldsymbol{\omega}$ is the vorticity vector

$$\boldsymbol{\omega} = \nabla \times \mathbf{v}, \qquad (32)$$

and $\mathbf{G}$ is the so-called hyper-stress. The associated generalized momentum balance reads:

$$\rho \dot{\mathbf{v}} = \nabla \cdot \mathbf{T} + \nabla \times (\nabla \cdot \mathbf{G}). \qquad (33)$$

The following constitutive equation was assumed in [25] for the fluid hyper-stress:

$$\mathbf{G} = \eta \mathcal{L}_f^2 (\nabla \boldsymbol{\omega} + \iota (\nabla \boldsymbol{\omega})^T), \qquad (34)$$

where the parameter $\iota \in [-1, 1]$ controls the asymmetry of the hyper-stress and ensures a non-negative dissipation. In the previous equation $\mathcal{L}_f^2 > 0$ is known as the fluid gradient length.

In the case of an incompressible fluid, the introduction of both the stress and the hyper-stress constitutive equations into the generalized momentum balance leads to:

$$\rho \dot{\mathbf{v}} = -\nabla p + \eta \Delta (\mathbf{v} - \mathcal{L}_f^2 \Delta \mathbf{v}). \qquad (35)$$

where $\eta \mathcal{L}_f^2$ represents the so-called hyper-viscosity [25].

See [25] for a discussion on the associated boundary conditions. The higher-order velocity derivatives involved in Eq. (35) require the enforcement of a larger number of boundary conditions compared to the standard first-gradient formulation.
<figure><img src="image_4.png" /><figcaption>Fig. 3 Hydrodynamic forces applied on a flexible rod immersed in a Newtonian fluid.</figcaption></figure>

## 4.2 Second-gradient description of a dilute suspension of flexible rods

Inspired by the use of the vorticity $\omega$ in the above second-gradient fluid flow description, we consider now the rod beads being subjected to two actions, as depicted in Fig. 3:

- A hydrodynamic force $\mathbf{F}$ that depends on the difference of velocities between the fluid and the bead. The fluid velocity at the bead position $\mathbf{v}(\mathbf{p}L)$ taking into account second-gradient effects reads

$$ \mathbf{v}(\mathbf{p}L) = \mathbf{v}_0 + \nabla \mathbf{v} \cdot \mathbf{p}L + \mathbf{H} : (\mathbf{p} \otimes \mathbf{p})L^2. \quad (36) $$

The bead velocity is again given by $\mathbf{v}_G + \dot{\mathbf{p}}L$. Thus, the resulting hydrodynamic force reads

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

- A hydrodynamic torque $\mathbf{M}$ that depends on the difference between the differential vorticity at the bead position and the bead rotary velocity $\dot{\varphi}$. The differential vorticity is the difference between the vorticity existing at the bead position minus the one existing at the rod center of gravity. This torque could have its origin in the distributed forces applied along the rod length as illustrated in Section 3. Thus, the resulting torque reads

$$ \mathbf{M}(\mathbf{p}L) = \xi^R (\nabla \omega \cdot \mathbf{p}L - \dot{\varphi}(\mathbf{p}L)), \quad (38) $$

where $\xi^R$ is the rotary friction coefficient.

A forces balance yields

$$ \mathbf{v}_0 - \mathbf{v}_G = -\mathbf{H} : (\mathbf{p} \otimes \mathbf{p})L^2, \quad (39) $$

implying a second-order relative velocity of the rod center of gravity with respect to the fluid velocity at this position. Again, for notational simplicity, we consider $\mathbf{F} = \mathbf{F}(\mathbf{p}L) = -\mathbf{F}(-pL)$.
As soon as we consider a first gradient of the vorticity-based bending mechanism $\nabla\omega \cdot \mathbf{p}L$, it is easy to prove (see Appendix C) that the rod kinematics remains unchanged relative to Jeffery's model, and that this term only affects rod bending. Thus, as proved in Appendix C, we obtain:

- The standard expression of the forces applying on the rod beads:

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

- Jeffery's rod kinematics:

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

- The bending mechanism. By considering the notation $\dot{\varphi} = \dot{\varphi}(\mathbf{p}L) = -\dot{\varphi}(-pL)$ and assuming a linear elastic behavior of the rod, the relation between the applied torque $\mathbf{M}$ and the bending angle $\varphi$ at the rod bead is given by

$$ \varphi = \frac{L}{EI} \mathbf{M}, \qquad (42) $$

or

$$ \mathbf{M} = \frac{EI}{L} \varphi = \kappa \varphi, \qquad (43) $$

where $E$ is the elastic modulus and $I$ the moment of area of the rod cross section. Introducing the expression of the moment $\mathbf{M}$ (38) into Eq. (43), we obtain

$$ \xi^R (\nabla\omega \cdot \mathbf{p}L - \dot{\varphi}) = \kappa \varphi, \qquad (44) $$

whose integration allows calculating $\varphi$ and from it the moment $\mathbf{M}$.

*Remark 6.* Tt is noteworthy to note that:

- In first-gradient flows, rods orient according to Jeffery's equation without experiencing bending effects because the vorticity gradient vanishes.

- In second-order flows, rods orient according to Jeffery's equation experiencing a bending induced by the second gradient.

- In rigid kinematics, rod moves with the fluid without experiencing any relative movement. Because the vorticity gradient vanishes there are not bending effects, ensuring the formulation objectivity.

### 4.2.1 Rod kinematics

The rod kinematics is then fully described by $\mathbf{p}$ and $\varphi$. Knowing $\mathbf{p}$, we can locate both rod beads. Then, knowing the beads bending angle $\varphi$ and all the forces and moments being applied at both ends, the deformed rod configuration is parabolic and symmetric with respect to the center of gravity. Thus, from $(\mathbf{p} \text{ and } \varphi)$, we can predict the bent configuration.

The evolution of both descriptors are given by Eqs. (41) and (44). For illustrating the behaviour, we consider Poiseuille's flow defined in Appendix A. A rod is tracked all along its pathline. In order to emphasize the representation, we make use of non-physical parameters in the orientation and bending models, Eqs. (41) and (44) respectively: $L = 1$, $\kappa = 4$ and $\xi^R = 1$.
<figure><img src="image_5.png" /><figcaption>Fig. 4 Evolution of a flexible rod immersed in a second-gradient Poiseuille flow.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Evolution of a flexible rod immersed in a second-gradient Poiseuille flow, as seen in the rod frame.</figcaption></figure>

The integration of Eqs. (41) and (44) yields $\mathbf{p}(t)$ and $\varphi(t)$, that allow us to depict the deformed rod at the position of its center of gravity $\mathbf{x}_G(t)$. Figure 4 shows the rod evolution. In order to emphasize the configuration evolution, we superpose in Fig. 5 all rod centers of gravity.

We can notice that, as dictated by the Jeffery equation, the rod rotates counter-clockwise. Moreover, the gradient of vorticity activates rod bending, mainly when that gradient is maximum $\theta = \frac{\pi}{2}$. Then the rod approaches the equilibrium steady-state orientation $\theta = \pi$, its rotary velocity decreases and bending relaxes because the vorticity gradient vanishes progressively when approaching the flow direction.

### 4.2.2 Induced stresses

The stress has two components, the standard one and the one related to the hyper-stress. The first one is related to forces applied at both opposite beads acting in the rod direction and includes the suspending medium contribution $\boldsymbol{\tau}^f = 2\eta\mathbf{D}$:

$$ \begin{gathered} \boldsymbol{\tau} = \boldsymbol{\tau}^f + \boldsymbol{\tau}^r = 2\eta\mathbf{D} + \sum_{i=1}^{N} \mathbf{p}_i \otimes \mathbf{F}_i = \\ 2\eta\mathbf{D} + \beta\nabla\mathbf{v} : \left( \sum_{i=1}^{N} \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \otimes \mathbf{p}_i \right), \end{gathered} \quad (45) $$

with the total stress $\mathbf{T}$ given by

$$ \mathbf{T} = -p \mathbf{I} + \boldsymbol{\tau}, \quad (46) $$
which, as expected, is symmetric.

The second one is related to second-gradient fluid contribution and rod bending. We consider a partition of the total hyper-stress **G** consisting of the fluid **G**^f and the rods **G**^r contributions, with **G** = **G**^f + **G**^r.

We consider a standard form of the second-gradient fluid contribution **G**^f,

$$ \mathbf{G}^f = \eta \mathcal{L}_f^2 (\nabla \omega + \iota (\nabla \omega)^T), \quad (47) $$

where $\mathcal{L}_f^2$ is the fluid gradient length.

On the other hand, we compute the rod contribution to the hyper-stress **G**^r from

$$ \mathbf{G}^r = \varrho \sum_{i=1}^{N} \mathbf{p}_i \otimes \mathbf{M}_i, \quad (48) $$

from which we assume the more general form

$$ \tilde{\mathbf{G}}^r = \mathbf{G}^r + \zeta (\mathbf{G}^r)^T. \quad (49) $$

## 4.3 Mesoscopic description

Making use of the orientation distribution function $\Psi(\mathbf{x}, t, \mathbf{p}, \varphi)$ (see Appendix D for a summary on micro-to-macro descriptions), we can proceed to the calculation of the stress within a continuous macroscopic description. When considering the stress expression (45) and the use of the quadratic closure relation $\mathbf{A} \approx \mathbf{a} \otimes \mathbf{a}$, we obtain:

$$ \mathbf{\tau} = \mathbf{\tau}^f + \mathbf{\tau}^r = 2\eta\mathbf{D} + 2\eta N_p (\mathbf{D} : \mathbf{A}) \approx 2\eta\mathbf{D} + 2\eta N_p (\mathbf{D} : \mathbf{a})\mathbf{a}, \quad (50) $$

with $\mathbf{T} = -\mathbf{p}\mathbf{I} + \mathbf{\tau}$.

On the other hand, the hyper-stress is given by:

$$ \mathbf{G}^r = \varrho \int_{S \times C} \mathbf{p} \otimes \mathbf{M} \Psi \, d\mathbf{p} \, d\varphi, \quad (51) $$

and, by using the expression of the moment,

$$ \mathbf{G}^r = \varrho\kappa \int_{S \times C} \mathbf{p} \otimes \varphi \Psi \, d\mathbf{p} \, d\varphi = \tilde{\kappa} \int_{S \times C} \mathbf{p} \otimes \varphi \Psi \, d\mathbf{p} \, d\varphi. \quad (52) $$

In the previous expressions, the distribution function $\Psi = \Psi(\mathbf{x}, t, \mathbf{p}, \varphi)$ (its dependence on the different coordinates is omitted for the sake of clarity) gives the fraction of rods that at position $\mathbf{x}$ and time $t$ have a conformation given by $(\mathbf{p}, \varphi)$. The domains $S$ and $C$ refer respectively to the domains in which coordinates $\mathbf{p}$ and $\varphi$ are defined.

*Remark 7.* One could expect that the moment (52) vanishes, however taking into account that $\Psi(\mathbf{p}, \varphi) = \Psi(-\mathbf{p}, -\varphi)$, the integral does not vanish.

In order to close the formulation, we consider the second-order tensor $\mathbf{g}$ (such that $\mathbf{G}^r = \tilde{\kappa} \mathbf{g}$):

$$ \mathbf{g} = \int_{S \times C} \mathbf{p} \otimes \varphi \Psi \, d\mathbf{p} \, d\varphi, \quad (53) $$
and compute its time derivative taking into account the expressions of $\dot{\mathbf{p}}$ (41) and $\dot{\varphi}$ (44):

$$ 
\begin{aligned}
\dot{\mathbf{g}} &= \int_{S \times C} (\dot{\mathbf{p}} \otimes \varphi + \mathbf{p} \otimes \dot{\varphi}) \Psi \, d\mathbf{p} \, d\varphi = \\
&\quad \int_{S \times C} (\nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}) \otimes \varphi \Psi \, d\mathbf{p} \, d\varphi + \\
&\quad \int_{S \times C} \mathbf{p} \otimes \left( \nabla \omega \cdot \mathbf{p} \mathcal{L} - \frac{\kappa}{\xi R} \varphi \right) \Psi \, d\mathbf{p} \, d\varphi = \\
&\quad \nabla \mathbf{v} \cdot \mathbf{g} - \nabla \mathbf{v} : \left( \int_{S \times C} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \varphi \Psi \, d\mathbf{p} \, d\varphi \right) + L\mathbf{a} \cdot (\nabla \omega)^T - \frac{\kappa}{\xi R} \mathbf{g}.
\end{aligned}
\quad (54) $$

If we adopt the following closure relation:

$$ \int_{S \times C} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \varphi \Psi \, d\mathbf{p} \, d\varphi \approx \mathbf{a} \otimes \mathbf{g}, \quad (55) $$

then we obtain a closed form for the evolution of g:

$$ \dot{\mathbf{g}} = \nabla \mathbf{v} \cdot \mathbf{g} - (\nabla \mathbf{v} : \mathbf{a})\mathbf{g} + L\mathbf{a} \cdot (\nabla \omega)^T - \frac{\kappa}{\xi R} \mathbf{g}. \quad (56) $$

Multiplying Eq. (56) by $\tilde{\kappa}$, we obtain the time evolution of the contribution of rods to the hyper-stress:

$$ \dot{\mathbf{G}}^r = \nabla \mathbf{v} \cdot \mathbf{G}^r - (\nabla \mathbf{v} : \mathbf{a})\mathbf{G}^r + L\tilde{\kappa}\mathbf{a} \cdot (\nabla \omega)^T - \frac{\kappa}{\xi R}\mathbf{G}^r. \quad (57) $$

If we consider the coefficient affecting the vorticity gradient in Eq. (57) and take into account the relations $\tilde{\kappa} = \rho\kappa$ and $\kappa = EI/L$, we obtain

$$ L\tilde{\kappa} = E(\rho I) = E\mathcal{L}_r^2, \quad (58) $$

where $\mathcal{L}_r^2$ represents the rod gradient length. The last term in equation (57) involves the coefficient $\kappa/\xi^R$ with reciprocal time units, which, in absence of flow, controls the relaxation to the undeformed reference configuration. Thus, the inverse of this coefficient has the meaning of a relaxation time that we denote by $\mathcal{T}$. Moreover, taking into account the symmetry of tensor $\mathbf{a}$, we have $\nabla\mathbf{v}:\mathbf{a} = \mathbf{D}:\mathbf{a}$. Thus, Eq. (57) can be rewritten as

$$ \dot{\mathbf{G}}^r = \nabla \mathbf{v} \cdot \mathbf{G}^r - (\mathbf{D} : \mathbf{a}) \mathbf{G}^r + E \mathcal{L}_r^2 \mathbf{a} \cdot (\nabla \omega)^T - \frac{1}{\mathcal{T}} \mathbf{G}^r. \quad (59) $$

It is important to notice that the above equation involves two closure relations, the one related to the fourth-order orientation tensor $\mathbf{A}$ and the one expressed in Eq. (55).

The macroscopic flow model can thus be summarized as follows:

- Generalized momentum balance:

$$ \rho \dot{\mathbf{v}} = \nabla \cdot \mathbf{T} + \nabla \times (\nabla \cdot \mathbf{G}), \quad (60) $$

- Mass balance

$$ \nabla \cdot \mathbf{v} = 0, \quad (61) $$
— Constitutive equation:

$$ T = -p I + \tau, \qquad (62) $$

with $\tau = \tau^f + \tau^r$:

$$ \tau^f = 2\eta D, \qquad (63) $$

and

$$ \tau^r = 2\eta N_p (\nabla \mathbf{v} : \mathbf{a}) \mathbf{a} = 2\eta N_p (\mathbf{D} : \mathbf{a}) \mathbf{a}, \qquad (64) $$

— Hyper-stress $\mathbf{G} = \mathbf{G}^f + \tilde{\mathbf{G}}^r$:

$$ \mathbf{G}^f = \eta \mathcal{L}_f^2 (\nabla \boldsymbol{\omega} + \iota (\nabla \boldsymbol{\omega})^T), \qquad (65) $$

and

$$ \tilde{\mathbf{G}}^r = \mathbf{G}^r + \zeta (\mathbf{G}^r)^T, \qquad (66) $$

with

$$ \dot{\mathbf{G}}^r = \nabla \mathbf{v} \cdot \mathbf{G}^r - (\mathbf{D} : \mathbf{a}) \mathbf{G}^r + E \mathcal{L}_r^2 \mathbf{a} \cdot (\nabla \boldsymbol{\omega})^T - \frac{1}{\tau} \mathbf{G}^r. \qquad (67) $$

## 5 Linear viscoelastic behavior predicted by the proposed theory

Complex fluids are tested in the linear viscoelastic regime by applying a small-amplitude oscillatory flow in order to evaluate the frequency dependence of the complex modulus. The latter has a real part, the so-called storage modulus $G'$, and an imaginary part called loss modulus, $G''$.

Obviously, in viscous fluids and elastic solids, the storage and loss modulus vanish respectively. In general, viscoelastic materials, both behaviours coexist.

In order to perform a linear viscoelastic analysis for standard (first-gradient) fluids, it suffices to apply a small-amplitude oscillatory simple shear flow, since only first-order derivatives of the velocity field are involved in the model. In the present second-gradient framework, however, that involves second-order derivatives, the small-amplitude oscillatory flow must contain a richer kinematics.

Thus, we enforce the following displacement field

$$ \delta = \begin{pmatrix} \delta_0 f(y) e^{i \omega t} \\ 0 \\ 0 \end{pmatrix}, \qquad (68) $$

where $f(y)$ is a function depending on the $y$-coordinate, that can be differentiated at least twice, and $\delta_0$ is the oscillation amplitude. The associate velocity field is given by

$$ \mathbf{v} = \begin{pmatrix} u \\ v \\ w \end{pmatrix} = \begin{pmatrix} i\omega\delta_0 f(y) e^{i\omega t} \\ 0 \\ 0 \end{pmatrix}, \qquad (69) $$

and its gradient reads

$$ \nabla \mathbf{v} = \begin{pmatrix} 0 & i\omega\delta_0 f'(y) e^{i\omega t} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (70) $$
where $f'(y) = \frac{df(y)}{dy}$.

The flow is applied on an initially isotropic suspension and we can ignore the flow-induced orientation (the small-amplitude deformation implies a negligible change of the orientation state). Thus, we can assume that the orientation tensor **a** remains isotropic at all times:

$$ \mathbf{a} = \frac{1}{3}\mathbf{I}, \qquad (71) $$

where **I** denotes the unit tensor.

The suspension constitutive equation involves a purely viscous contribution,

$$ \tau = \tau^f + \tau^r = 2\eta\mathbf{D} + 2\eta N_p(\mathbf{D} : \mathbf{a})\mathbf{a}, \qquad (72) $$

that vanishes as soon as the flow stops and that consequently only concerns the loss modulus $G''$ (i.e. it is not involved in the elastic behavior).

For this reason, in what follows, we focus on the second contribution to the stress, more precisely on the hyper-stress **G**, $\mathbf{G} = \mathbf{G}^f + \tilde{\mathbf{G}}^r$, particularly on the one associated to the presence of rods $\tilde{\mathbf{G}}^r$, because the one related to the fluid $\mathbf{G}^f$

$$ \mathbf{G}^f = \eta \mathcal{L}_f^2 (\nabla \boldsymbol{\omega} + \iota (\nabla \boldsymbol{\omega})^T), \qquad (73) $$

is purely viscous, thus vanishing as soon as the flow stops.

When focusing on the rod contribution to the hyper-stress

$$ \tilde{\mathbf{G}}^r = \mathbf{G}^r + \zeta (\mathbf{G}^r)^T, \qquad (74) $$

with

$$ \dot{\mathbf{G}}^r = \nabla \mathbf{v} \cdot \mathbf{G}^r - (\mathbf{D} : \mathbf{a}) \mathbf{G}^r + E \mathcal{L}_r^2 \mathbf{a} \cdot (\nabla \boldsymbol{\omega})^T - \frac{1}{T} \mathbf{G}^r, \qquad (75) $$

we can neglect the first term of the right-hand side, because it only ensures objectivity of the evolution equation in the non-linear regime. Thus, the linearized problem reads

$$ \dot{\mathbf{G}}^{r,lin} = -(\mathbf{D} : \mathbf{a}) \mathbf{G}^{r,lin} + E \mathcal{L}_r^2 \mathbf{a} \cdot (\nabla \boldsymbol{\omega})^T - \frac{1}{T} \mathbf{G}^{r,lin}. \qquad (76) $$

In Eq. (76) the term $(\mathbf{D} : \mathbf{a}) = 0$ and the vorticity $\boldsymbol{\omega}$ is given by

$$ \boldsymbol{\omega} = \begin{pmatrix} 0 \\ 0 \\ -i\omega\delta_0 f'(y)e^{i\omega t} \end{pmatrix}, \qquad (77) $$

which implies

$$ \nabla \boldsymbol{\omega} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & -i\omega\delta_0 f''(y)e^{i\omega t} & 0 \end{pmatrix}, \qquad (78) $$

and consequently,

$$ \mathbf{a} \cdot (\nabla \boldsymbol{\omega})^T = \begin{pmatrix} 0 & 0 & 0 \\ 0 & -\frac{1}{3}i\omega\delta_0 f''(y)e^{i\omega t} & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & -\frac{1}{3}i\omega\delta_0 f''(y) & 0 \\ 0 & 0 & 0 \end{pmatrix} e^{i\omega t} = \mathbf{A}e^{i\omega t}. \qquad (79) $$
We assume that the response $G^{r,lin}$ to the enforced oscillation has the same frequency (in view of the linearity of the model when assuming small-amplitude oscillations) but can exhibit a phase delay $\varphi$, that is

$$G^{r,lin} = G_0 e^{i\omega t - i\varphi} = G_0 e^{i\omega t}, \qquad (80)$$

where the complex amplitude $G_0$ is such that $G_0 = G_0 e^{-i\varphi}$.

Introducing expressions (79) and (80) into Eq. (76), we obtain

$$G_0 i\omega e^{i\omega t} = \alpha A e^{i\omega t} - \frac{1}{\mathcal{T}} G_0 e^{i\omega t}, \qquad (81)$$

or

$$G_0 \left( i\omega + \frac{1}{\mathcal{T}} \right) = \alpha A, \qquad (82)$$

with $\alpha = E\mathcal{L}_r^2$.

Algebraic manipulations in Eq. (82) yield

$$G_0 = \frac{\mathcal{T}\alpha(1 - i\omega\mathcal{T})A}{1 + \omega^2\mathcal{T}^2}. \qquad (83)$$

The only non-zero component of $G_0$ is $G_{0,23}$, which by taking into account the expression of $A_{23}$ from Eq. (79) reads

$$G_{0,23} = \frac{\beta \mathcal{T} \omega^2 + i\beta\omega}{1 + \omega^2 \mathcal{T}^2}, \qquad (84)$$

with $\beta = -\frac{\alpha \mathcal{T} \delta_0 f''(y)}{3}$. When considering a parabolic viscoemtric flow, we have $f''(y) = \text{cst.}$

*Remark 8.* The contribution of the hyper-stress to the tension on a plate surface with unit normal $n$ is given by $n \times \nabla \cdot \tilde{G}^r$. For a surface defined by $n^T = (0, 1, 0)$, the only component that resists the flow $v^T = (u, 0, 0)$ is $\tilde{G}_{32}^r$, that results from $G_{23}^r$ and $G_{32}^r$, the first of them being in the present case non-zero.

We can conclude from the above analysis that the storage modulus $G' = \frac{\Re(G_{0,23})}{\delta_0}$, responsible for the elastic behaviour, has a dependence on the frequency $\omega$, exhibiting a slope of 2 at low frequencies and becoming constant (plateau) at high frequencies, i.e. the well-known Maxwellian behaviour observed in many complex fluids.

# 6 Inception of shear flow followed by relaxation

In this section, we wish to visualize from a physical point of view the subtle coupling between flow and rod elasticity, by considering two scenarios: (i) first, a prescribed shear flow acting on initially-straight flexible rods and yielding rod bending, followed by (ii) the progressive cessation of flow once the bent rods relax to their straight configurations after the flow source (pressure drop) has been abruptly set to zero.
## 6.1 Inception of shear flow

We prescribe the shear flow

$$ \mathbf{v} = \begin{pmatrix} u \\ v \\ w \end{pmatrix} = \begin{pmatrix} f(y) \\ 0 \\ 0 \end{pmatrix}, \qquad (85) $$

with $f(y)$ at least four times differentiable since the second-gradient flow model involves the curl of the hyper-stress divergence, $\nabla \times \nabla \cdot \mathbf{G}$, with the hyper-stress depending on the vorticity gradient (see Eqs. (33) and (34)).

We focus on a rod located on the streamline $y_0 = \text{cst}$, such that $f'(y_0) = \frac{du}{dy}|_{y_0} < 0$.

The velocity gradient, vorticity and vorticity gradient read respectively

$$ \nabla \mathbf{v} = \begin{pmatrix} 0 & f'(y_0) & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (86) $$

$$ \boldsymbol{\omega} = \begin{pmatrix} 0 \\ 0 \\ -f'(y_0) \end{pmatrix}, \qquad (87) $$

and

$$ \nabla \boldsymbol{\omega} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & -f''(y_0) & 0 \end{pmatrix}. \qquad (88) $$

We assume that at the initial time $t = 0$ the rods are fully contained in the $x-y$ plane. Thus, the orientation tensor has the form

$$ \mathbf{a}(t=0) = \mathbf{a}_0 = \begin{pmatrix} \mathbf{a}_{0,11} & \mathbf{a}_{0,12} & 0 \\ \mathbf{a}_{0,21} & \mathbf{a}_{0,22} & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (89) $$

where $\mathbf{a}_{0,12} = \mathbf{a}_{0,21}$ and $\mathbf{a}_{0,22} = 1 - \mathbf{a}_{0,11}$.

The hyper-stress evolution equation reads

$$ \dot{\mathbf{G}}^r = \nabla \mathbf{v} \cdot \mathbf{G}^r - (\mathbf{D} : \mathbf{a}) \mathbf{G}^r + E \mathcal{L}_r^2 \mathbf{a} \cdot (\nabla \boldsymbol{\omega})^T - \frac{1}{T} \mathbf{G}^r. \qquad (90) $$

At $t = 0$, assuming a fully relaxed system (straight rods), $\mathbf{G}^r(t = 0) = \mathbf{G}_0^r = \mathbf{0}$, the hyper-stress evolution only depends on the term $E \mathcal{L}_r^2 \mathbf{a}_0 \cdot (\nabla \boldsymbol{\omega})^T$, i.e.

$$ \dot{\mathbf{G}}^r(t=0) \propto \mathbf{a}_0 \cdot (\nabla \boldsymbol{\omega})^T, \qquad (91) $$

which, taking into account Eqs. (88) and (89), gives

$$ \dot{\mathbf{G}}_0^r \propto \begin{pmatrix} 0 & 0 & -\mathbf{a}_{0,12} f''(y_0) \\ 0 & 0 & -\mathbf{a}_{0,22} f''(y_0) \\ 0 & 0 & 0 \end{pmatrix}. \qquad (92) $$

We consider two scenarios:
- All rods are initially aligned in the flow direction, implying

$$ \mathbf{a}_0 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (93) $$

and according to Eq. (92) $\dot{\mathbf{G}}_0^r = \mathbf{0}$. In this case, rods remain oriented in the flow direction (according to Jeffery's equation) and they do not experience bending because $\mathbf{G}^r(t) = \mathbf{G}_0^r = \mathbf{0}$.

- All rods are initially aligned perpendicularly to the flow direction, implying

$$ \mathbf{a}_0 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (94) $$

and thus according to Eq. (92),

$$ \dot{\mathbf{G}}_0^r \propto \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -f''(y_0) \\ 0 & 0 & 0 \end{pmatrix}. \qquad (95) $$

Consequently, the hyper-stress is given at time $\delta t$ by

$$ \mathbf{G}^r(t = \delta t) = \mathbf{G}_\delta^r \approx E \mathcal{L}_r^2 \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -f''(y_0)\delta t \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & \mathbf{G}_{\delta,23}^r \\ 0 & 0 & 0 \end{pmatrix}. \qquad (96) $$

In order to interpret the physics of this result, we come back to the microscopic definition of the hyper-stress related to the rods (Eq. (52)), all of them being aligned initially in the direction given by the unit vector $\mathbf{p}_0$

$$ \mathbf{p}_0 = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}, \qquad (97) $$

which implies

$$ \tilde{\kappa} \mathbf{p}_0 \otimes \varphi = \mathbf{G}_\delta^r, \qquad (98) $$

and thus

$$ \varphi = \begin{pmatrix} 0 \\ 0 \\ \varphi \end{pmatrix}, \qquad (99) $$

with $\tilde{\kappa}\varphi = \mathbf{G}_{\delta,23}^r$. With $f''(y_0) < 0$, we have $\mathbf{G}_{\delta,23}^r > 0$ and then $\varphi > 0$ inducing rod bending as illustrated in Fig. 4.

We can go one step forward, integrating Eq. (92) from $t = \delta t$ to $t = 2\delta t$.

Now, $\mathbf{G}_\delta^r$ is given by Eq. (96). Integration of Jeffery's equation results in an updated orientation $\mathbf{p}(t = \delta t) = \mathbf{p}_\delta$ given by

$$ \mathbf{p}_\delta = \begin{pmatrix} \mathbf{p}_{\delta,x} \\ \mathbf{p}_{\delta,y} \\ 0 \end{pmatrix}, \qquad (100) $$

where $\mathbf{p}_{\delta,y} \approx 1$ and $\mathbf{p}_{\delta,x} < 0$, with $||\mathbf{p}_\delta|| = 1$.

Now the term $\mathbf{a}_\delta \cdot (\nabla\omega)^T$ activates the evolution of $\mathbf{G}_{13}^r$ and $\mathbf{G}_{23}^r$ and the objectivity of the resulting evolution is ensured by the term $\nabla\mathbf{v} \cdot \mathbf{G}_\delta^r$ in Eq. (92) that takes into account that the rod is rotating counterclockwise, in order not to include the rod rotation in the bead bending angle.
## 6.2 Progressive cessation of flow induced by rod relaxation

We consider now the state reached by the system at time $\delta t$ with all rods initially aligned perpendicularly to the flow direction (scenario discussed in the second item of the previous section).

In that case the hyper-stress reads

$$ \mathbf{G}_{\delta}^{r} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & \mathbf{G}_{\delta,23}^{r} \\ 0 & 0 & 0 \end{pmatrix}. \qquad (101) $$

Neglecting $\mathbf{G}^f$ (i.e. putting $\mathcal{L}_f = 0$) and for the sake of clarity neglecting also the contribution of rods to the stress (i.e. considering $N_p \approx 0$), the flow model (38) reduces to

$$ \frac{\partial p}{\partial x} = \eta \frac{\partial^2 u}{\partial y^2} + \zeta \frac{\partial^2 \mathbf{G}_{23}^r}{\partial y^2}. \qquad (102) $$

Thus, if we remove suddenly the external action applied for creating the flow, ie. the pressure drop $\frac{\partial p}{\partial x}$ is set to zero, the flow persists until the full relaxation of the hyper-stress, that is until rods recover their initially-straight configuration. In absence of viscosity, this relaxation occurs instantaneously, but the presence of the viscous fluid requires that flow occurs for accommodating relaxation of rod bending.

This mechanism, that can be physically fully interpreted, is due to the transient nature of the equation that governs the rod contribution to the hyper-stress. Assuming $\mathbf{G}^f \neq 0$ and $\mathbf{G}^r = 0$, as soon as the external action creating the flow is removed, the flow stops suddenly, because the fluid hyper-stress is not concerned by a relaxation mechanism, that is, it relaxes instantaneously, like a viscous stress.

## 7 Complex flow simulation

In this section, we solve the proposed suspension flow model in the driven cavity flow problem. The suspension is confined in a cavity $\Omega = [0, L]^2$, with $L = 2$. The velocity is specified at the upper wall as $\mathbf{v}^T(x, y = L) = (Lx - x^2)$, and it is set to zero in the remaining part of the domain boundary (left, right and bottom walls). Moreover, as proposed in [25], we enforce the extra-condition $\mathbf{n} \times \mathbf{G} \cdot \mathbf{n} = -e\tilde{t}al\omega \times \mathbf{n}$ in the whole domain boundary $\partial\Omega$, where $l$ and $\tilde{l}$ are respectively the so-called adherence length and the boundary viscosity.

In order to take the incompressibility constraint $\nabla \cdot \mathbf{v} = 0$ into account, we make use of the stream function formulation. The stream function $\Xi$ is related to the components of the velocity vector according to

$$ \left\{ \begin{array}{l} u = -\frac{\partial \Xi}{\partial y} \\ v = \frac{\partial \Xi}{\partial x} \\ w = 0 \end{array} \right. . \qquad (103) $$

Both velocity and vorticity are written in terms of the stream function $\Xi$. It is easy to verify that

$$ \omega = \begin{pmatrix} 0 \\ 0 \\ \Delta\Xi \end{pmatrix}. \qquad (104) $$
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Second-gradient model solution: flow streamlines.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Second-gradient model solution: rod bending angle at two different times $t_1$ (left) and $t_2$ (right), with $t_2 > t_1$.</figcaption></figure>

Due to the high order of the resulting partial differential equation, a spectral Chebyshev collocation technique is used here for discretizing the motion equation.

The computed streamlines related to the second-gradient flow model solution are depicted in Fig. 6. In this simulation, we considered $\mathcal{L}_f = 0$, $N_p \approx 0$ and $\mathcal{L}_r \ll 1$, which implies that the flow kinematics are very close to the ones related to a Newtonian fluid. Our main interest when considering the second-gradient description concerns the evaluation of the rod elastic loading, more than the analysis of the perturbed flow kinematics analyzed in detail in [25].

Figure 7 shows the absolute value of the rod bending angle $\varphi$ at two different times. The initial orientation state is full alignment in the y-direction everywhere in the computational domain $\Omega$.
## 8 Conclusions

We have developed first and second-gradient models for the kinematics of rods immersed in a Newtonian fluid, with the rods being modeled as dumbbells. It has been shown that in both descriptions the rod kinematics are described by Jeffery’s equation and that, with forces acting along the rod direction, rod bending cannot be activated. We have then considered an alternative approach wherein hydrodynamic forces are distributed all along the rod length. In that case, the rod kinematics are again described by Jeffery’s equation but, as soon as a second-gradient description is considered, rod bending is activated. From this microscopic description, we have derived a second-gradient macroscopic model that couples the rod behaviour with the flow kinematics. The proposed theory has been shown to predict a Maxwell linear viscoelastic behaviour. We have studied the model predictions in the inception and cessation of simple shear flow. Finally, we have illustrated the use of the proposed model in the simulation of complex flow in a driven cavity, assuming coupling between rod and flow kinematics.

A deeper analysis of the flow-induced bending and the effects of rod bending on flow kinematics is work in progress. The proposed model is also being extended to concentrated suspensions involving entangled flexible rods wherein the rotary drag is created by the rotation of the rod network. This extension of the proposed theory could constitute a valuable approach for modeling complex flows of concentrated suspensions involving flexible rods as encountered in SMC processes in composite manufacturing.

## A On the relative velocity between the fluid and the rod center of gravity

We showed in Section 2.2 that, when considering a second-gradient description, the rod center of gravity has a velocity different from the local fluid velocity. Even though one could think that this relative velocity implies rod migration, we here show that it is not the case. The relative velocity derived in Section 2.2 is given by

$$v_0 - v_G = -(\mathbf{H} : (\mathbf{p} \otimes \mathbf{p}))L^2. \qquad (105)$$

Let us consider a 2D Poiseuille flow described by the following kinematics:

$$v = \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \dot{\gamma} (H^2 - y^2) \\ 0 \end{pmatrix}, \qquad (106)$$

with $y \in [-H, H]$.

In this case, the only non-zero component of tensor $\mathbf{H}$ is:

$$\mathbf{H}_{122} = \frac{1}{2} \frac{\partial^2 u}{\partial^2 y} = -\dot{\gamma}, \qquad (107)$$

which acts on $(\mathbf{p} \otimes \mathbf{p})_{22}$. With $\mathbf{p}^T = (\cos \theta, \sin \theta)$, we obtain

$$v_0 - v_G = \begin{pmatrix} \dot{\gamma}L^2 \sin^2 \theta \\ 0 \end{pmatrix}. \qquad (108)$$

This shows that rods do not move towards regions of lower shear rates. They never leave their trajectories, $y = \text{cst}$, but depending on their orientation, their center of gravity moves faster or slower than the fluid.
$$ \mathbf{F} = \xi (\mathbf{v}_0 + \nabla \mathbf{v} \cdot \mathbf{L}\mathbf{p} + \mathbf{H} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{L}^2 + \mathbf{J} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}))\mathbf{L}^3 - \mathbf{L}\dot{\mathbf{p}}) $$

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Hydrodynamic forces acting on a rod immersed in a Newtonian fluid, considering the third-order velocity gradient.</figcaption></figure>

## B Third-gradient description

We assume here that the forces **F** that act on each bead depend again on the difference of velocities between the fluid and the bead, but the former now including the third-order velocity gradient **J**, i.e. that implies the fluid velocity at the bead position given by **v**(pL) = **v**₀ + ∇**v** · pL + (H : (p ⊗ p))L² + (J : (p ⊗ p ⊗ p))L³ and the second one by **v**<sub>G</sub> + **p**·L. Thus, the force **F**(pL) reads (see Fig. 8):

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

where the third-order velocity gradient **J** is the fourth-rank tensor with components **J**<sub>ijkm</sub> = 1/6 ∂³**v**<sub>i</sub> / ∂x<sub>j</sub>∂x<sub>k</sub>∂x<sub>m</sub>.

Obviously, if **F** acts on the bead pL, then for the opposite bead at −pL the resulting force reads

$$ \mathbf{F}(-\mathbf{p}L) = \xi (\mathbf{v}_0 - \nabla \mathbf{v} \cdot \mathbf{p}L + (\mathbf{H} : (\mathbf{p} \otimes \mathbf{p}))L^2 - (\mathbf{J} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}))L^3 - \mathbf{v}_G + \dot{\mathbf{p}}L). \quad (110) $$

By adding Eqs. (109) and (110) and enforcing the force balance

$$ \mathbf{F}(\mathbf{p}L) + \mathbf{F}(-\mathbf{p}L) = \mathbf{0}, \quad (111) $$

we obtain again **v**₀ − **v**<sub>G</sub> = −(H : (p ⊗ p))L².

As the resulting torque must also vanish, the only possibility is that the force **F** acts along **p**, that is **F** = λ**p**, with λ ∈ ℝ. Thus, we can write

$$ \lambda \mathbf{p} = \xi (\nabla \mathbf{v} \cdot \mathbf{p}L + (\mathbf{J} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}))L^3 - \dot{\mathbf{p}}L), \quad (112) $$

which multiplyed by **p** yields an expression of λ:

$$ \lambda = \xi L (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}) + L^2 (\mathbf{J} :: (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}))) \quad (113) $$

We thus have

$$ \mathbf{F} = \lambda \mathbf{p} = \xi L (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})\mathbf{p} + L^2 (\mathbf{J} :: (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p})))\mathbf{p}), \quad (114) $$

which yield the rotary velocity

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p} + L^2 (\mathbf{J} : (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p})) - (\mathbf{J} :: (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p})))\mathbf{p}, \quad (115) $$
where we can identify Jeffery's contribution $\dot{\mathbf{p}}^J$ and a third-order correction $\dot{\mathbf{p}}^H$ that is affected by $L^2$:

$$ \dot{\mathbf{p}}^J = \nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})) \mathbf{p}, \qquad (116) $$

and

$$ \dot{\mathbf{p}}^H = \mathbf{J} \cdot (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}) - (\mathbf{J} \cdot (\mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}))\mathbf{p}, \qquad (117) $$

with the total velocity given by:

$$ \dot{\mathbf{p}} = \dot{\mathbf{p}}^J + L^2 \dot{\mathbf{p}}^H. \qquad (118) $$

## C Kinematics of rods subjected to linear and angular drags

By enforcing the moment balance, we obtain

$$ 2L\mathbf{p} \times \mathbf{F} + \mathbf{M}(\mathbf{p}L) + \mathbf{M}(-\mathbf{p}L) = \mathbf{0}, \qquad (119) $$

or

$$ 2\xi L^2 (\mathbf{p} \times \nabla \mathbf{v} \cdot \mathbf{p} - \mathbf{p} \times \dot{\mathbf{p}}) - \xi^R (\dot{\varphi}(\mathbf{p}L) + \dot{\varphi}(-\mathbf{p}L)) = \mathbf{0}. \qquad (120) $$

We define $\mathbf{W}$ such that

$$ \dot{\mathbf{p}} = \mathbf{W} \times \mathbf{p}, \qquad (121) $$

and since $\mathbf{p}$ and $\mathbf{W}$ are perpendicular (for the sake of simplicity, we consider rods with a planar orientation), we have $\mathbf{p} \times \dot{\mathbf{p}} = \mathbf{W}$. Moreover, if we define the average angular velocity $\dot{\bar{\varphi}}$ as

$$ \dot{\bar{\varphi}} = \frac{\dot{\varphi}(\mathbf{p}L) + \dot{\varphi}(-\mathbf{p}L)}{2}, \qquad (122) $$

then Eq. (120) can be rewritten as:

$$ 2\xi L^2 (\mathbf{p} \times \nabla \mathbf{v} \cdot \mathbf{p} - \mathbf{W}) - 2\xi^R \dot{\bar{\varphi}} = \mathbf{0}, \qquad (123) $$

from which we obtain the rod rotary velocity

$$ \mathbf{W} = \mathbf{p} \times \nabla \mathbf{v} \cdot \mathbf{p} - \frac{\xi^R}{\xi L^2} \dot{\bar{\varphi}}. \qquad (124) $$

Taking into account Eq. (121) and the triple vector product property $\mathbf{a} \times \mathbf{b} \times \mathbf{c} = \mathbf{b}(\mathbf{a} \cdot \mathbf{c}) - \mathbf{c}(\mathbf{a} \cdot \mathbf{b})$, we have

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p}))\mathbf{p} - \frac{\xi^R}{\xi L^2} \dot{\bar{\varphi}} \times \mathbf{p}. \qquad (125) $$

Now, coming back to Eq. (37) and taking Eq. (125) into account,

$$ \mathbf{F} = \xi L \nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})\mathbf{p} + \frac{\xi^R}{L} \dot{\bar{\varphi}} \times \mathbf{p}. \qquad (126) $$

As we can notice, the force acting on the beads has a component in the rod direction given by

$$ \mathbf{F}^{\|} = \xi L \nabla \mathbf{v} : (\mathbf{p} \otimes \mathbf{p})\mathbf{p}, \qquad (127) $$

and another component perpendicular to it and contained in the plane where the rod orientation is defined:

$$ \mathbf{F}^{\perp} = \frac{\xi^R}{L} \dot{\bar{\varphi}} \times \mathbf{p}. \qquad (128) $$

Since $\mathbf{F}^{\perp}(\mathbf{p}L) = -\mathbf{F}^{\perp}(-\mathbf{p}L)$ and making use of symmetry considerations, equality of torques and bending angles at both rod beads are expected, that is:

$$ \mathbf{M}(\mathbf{p}L) = \mathbf{M}(-\mathbf{p}L), \qquad (129) $$

and

$$ \varphi(\mathbf{p}L) = \varphi(-\mathbf{p}L), \qquad (130) $$
respectively. Thus, we find $\bar{\phi} = \dot{\phi}$.

This result, however, is not compatible with

$$ \left\{ \begin{aligned} \mathbf{M}(\mathbf{p}L) &= \xi^R(\nabla\omega \cdot \mathbf{p}L - \dot{\phi}(\mathbf{p}L)) \\ \mathbf{M}(-\mathbf{p}L) &= \xi^R(-\nabla\omega \cdot \mathbf{p}L - \dot{\phi}(\mathbf{p}L)) \end{aligned} \right. . \qquad (131) $$

Thus, the only possibility is that $\dot{\phi}(\mathbf{p}L) = -\dot{\phi}(-\mathbf{p}L)$, implying $\mathbf{F}^{\perp}(\mathbf{p}L) = \mathbf{F}^{\perp}(-\mathbf{p}L) = \mathbf{0}$ (which yields the standard Jeffery solution for the rod rotary velocity $\dot{\mathbf{p}}$) and $\mathbf{M}(\mathbf{p}L) = -\mathbf{M}(-\mathbf{p}L)$, which fully agrees with Eq. (131).

# D From mesoscopic to macroscopic descriptions

## D.1 Kinetic theory description

In view of the very large number of rods involved in a typical suspension, the description based on the tracking of each individual particle, 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 rod individuality is lost in favour of a distribution function.

For the sake of simplicity, we consider rigid rods in what follows. Thus, one could describe the microstructure at a certain point $\mathbf{x}$ and time $t$ from the orientation distribution function $\Psi(\mathbf{x}, t, \mathbf{p})$ that gives 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 \Psi(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p} = 1, \quad \forall \mathbf{x}, \forall t, \qquad (132) $$

where $S$ 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 possible rod orientations.

Thus, we have substituted the fine-grain microscopic description, that requires the computation of each individual rod orientation, with the scalar multidimensional function $\Psi$. However, in order to use it, one needs to derive the equation governing the evolution of the orientation distribution function $\Psi$.

The balance ensuring conservation of orientation probability reads:

$$ \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\dot{\mathbf{x}} \Psi) + \frac{\partial}{\partial \mathbf{p}} (\dot{\mathbf{p}} \Psi) = 0, \qquad (133) $$

while, for inertialess rods and considering a first-gradient description $\dot{\mathbf{x}} = \mathbf{v}(\mathbf{x}, t)$, the rod rotary velocity is given by Jeffery's equation:

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

Eq. (133), known as the 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 description for the flow-induced orientation of individual rods given by the Jeffery equation (134).

The price to pay is the increase in the model dimensionality, since the orientation distribution is defined in a high-dimensional domain of dimension 5 in the general 3D case, i.e. $\Psi : (\mathbf{x}, t, \mathbf{p}) \to \mathbb{R}^+$ where $\mathbf{x} \in \Omega \subset \mathbb{R}^3, t \in \mathcal{I} \subset \mathbb{R}^+, \mathbf{p} \in S$.

The extra-stress due to the presence of such a population of rods is determined by adding the contribution of each rod to the stress:

$$ \begin{aligned} \tau &= \tau^f + \tau^r = \\ & 2\eta \mathbf{D} + 2\eta N_p \nabla \mathbf{v} : \left( \int_S \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \Psi(\mathbf{x}, t, \mathbf{p}) \, d\mathbf{p} \right) = \\ & 2\eta \mathbf{D} + 2\eta N_p (\nabla \mathbf{v} : \mathbf{A}), \end{aligned} \qquad (135) $$
where **A** represents the fourth-order moment of the orientation distribution function

$$ \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 (136) $$

and $N_p$ is the so-called particle number that depends on rod concentration and on the friction coefficient. In view of the symmetry of **A**, the gradient of velocities in the previous expression can be replaced by the rate of strain tensor **D**:

$$ \tau = \tau^f + \tau^r = 2\eta\mathbf{D} + 2\eta N_p (\mathbf{D} : \mathbf{A}). \qquad (137) $$

## D.2 Macroscopic description

Fokker-Planck based descriptions are rarely used in practice in view of the curse of dimensionality that the introduction of conformation coordinates implies (here, the rod orientation). Thus, standard mesh-based discretization techniques, such as finite differences, finite elements or finite volumes, fail when addressing models defined in high-dimensional spaces. For this reason, mesoscopic models are usually coarsened one step further in order to derive macroscopic models defined in standard physical domains, involving only space and time coordinates.

In this section, we illustrate the transition from the mesoscopic to the macroscopic scale. At the macroscopic scale, moments of the orientation distribution are used for describing the microstructure.

We consider the second moment **a**:

$$ \mathbf{a} = \int_S \mathbf{p} \otimes \mathbf{p} \Psi \, d\mathbf{p} \qquad (138) $$

whose time derivative reads:

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

A closure relation is needed in order to express the fourth-order moment **A** as a function of the lower-order moments (as odd moments vanish because of the symmetry of $\Psi$, the only non-zero lower moment is the second-order moment **a**). Different closure relations have been introduced and widely used [2] [13] [26] [31]. For example, in the quadratic closure relation (that is exact only when all rods are locally aligned in the same direction), the fourth-order moment is approximated as

$$ \mathbf{A} \approx \mathbf{a} \otimes \mathbf{a}. \qquad (140) $$

We thus have

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

and the stress tensor is given by

$$ \tau = \tau^f + \tau^r = 2\eta\mathbf{D} + 2\eta N_p (\mathbf{D} : \mathbf{A}) \approx 2\eta\mathbf{D} + 2\eta N_p (\mathbf{D} : \mathbf{a}) \mathbf{a}. \qquad (142) $$
References

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

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

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

4. J. Azaiez, K. Chiba, F. Chinesta, A. Poitou. State-of-the-art on numerical simulation of fiber-reinforced thermoplastic forming processes. Archives of Computational Methods in Engineering, **9/2**, 141-198, (2002).

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

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

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

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

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

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

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

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

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

14. C. Eringen, Nonlocal continuum field theories, Springer 2002.

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

16. E. Fried, M.E. Gurtin, Traction, balances, and boundary conditions for non-simple materials with application to liquid flow at small length scales, Arch. Ration. Mech. Anal., **182**, 513-554, 2006.

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

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**, 683-712, 1972.

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

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

21. J. Hinch, The distortion of a flexible inextensible thread in a shearing flow, Fluid Mech., **74**, 317-333, 1976.

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

23. R. Keunings. On the Peterlin approximation for finitely extensible dumbells, J. Non-Newtonian Fluid Mech., **68**, 85-100, 1997.

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

25. T. Kim a, J. Dolbow, E. Fried. A numerical method for a second-gradient theory of incompressible fluid flow, Journal of Computational Physics **223**, 551-570, 2007.

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

27. 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.
28. A. Ma, F. Chinesta, A. Ammar, M. Mackley, Rheological modelling of Carbon Nanotube aggregate suspensions, Journal of Rheology, **52/6**, 1311-1330, 2008.

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

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

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

32. R. Shanker, J.W. Gillespie, S.I. Güçeri, On the effect of nonhomogeneous flow fields on the orientation distribution and rheology of fiber suspensions, Polymer Engineering and Science, **31**, 161-171, 1991.

33. R. Shanker, The effect of non homogeneous flow fields and hydrodynamic interactions on the rheology of fiber suspensions, PhD dissertation, University of Delaware, 1991.

34. U. Strautins, A. Latz, Flow-driven orientation dynamics of semiflexible fiber systems, Rheol Acta, **46**, 1057-1064, 2007.

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