# A numerical analysis of the evolution of bundle orientation in concentrated fibre-bundle suspensions

Pierre J.J. Dumont, Steven Le Corre, Laurent Orgéas, Denis Favier

► To cite this version:

Pierre J.J. Dumont, Steven Le Corre, Laurent Orgéas, Denis Favier. A numerical analysis of the evolution of bundle orientation in concentrated fibre-bundle suspensions. Journal of Non-Newtonian Fluid Mechanics, 2009, 160 (2-3), pp.76-92. 10.1016/j.jnnfm.2009.03.001. hal-01004927

HAL Id: hal-01004927

https://hal.science/hal-01004927

Submitted on 11 Jun 2014

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 numerical analysis of the evolution of bundle orientation in concentrated fibre-bundle suspensions

P.J.J. Dumont<sup>a</sup>, S. Le Corre<sup>b,*</sup>, L. Orgéas<sup>c</sup>, D. Favier<sup>c</sup>

<sup>a</sup> CNRS/Universités de Grenoble, Laboratoire de Génie des Procédés Papetiers (LGP2), 461 rue de la Papeterie, BP 65, 38402 Saint-Martin-d'Hères cedex, France

<sup>b</sup> GeM - Institut de Recherche en Génie Civil et Mécanique, CNRS/École Centrale de Nantes, BP 92101, 44321 Nantes cedex 3, France

<sup>c</sup> CNRS/Universités de Grenoble, Laboratoire Sols-Solides-Structures-Risques (3S-R), BP 53, 38041 Grenoble cedex, France

The rheology of Newtonian concentrated fibre-bundle suspensions with nearly planar orientation states is investigated within the framework of the homogenization method for periodic discrete structures. These suspensions are seen as forming a connected network of bundles. At their contact points, the flow induces local Newtonian interaction forces and moments. "Numerical rheometry experiments" are performed on representative elementary volumes of these suspensions to explore the influence of bundle volume fraction and orientation state on the rate of change of bundle orientation. Three bundle orientation distributions are particularly investigated: Gaussian-based, uniform and crenellated orientation distributions. Two types of numerical simulations are performed, i.e., instantaneous and time-evolution ones. The obtained numerical results permit to discuss the role of the shape of the orientation distribution, typically on the rate of change of the second-order orientation tensor. They are also compared to the well-known Jeffery theory and its subsequent modifications like the Folgar-Tucker theory. The accuracy of some well-known closure approximation functions of the fourth-order orientation tensor is examined. It appears that, in the case treated here where the suspending fluid and the fibre-fibre interactions are Newtonian, the Jeffery's equation gives a fairly good fit (the best fit of all examined theories) of the homogenization results, despite the completely different physics of both approaches. Some differences are observed in situations where fibre-bundles are highly aligned along the flow direction. It is also observed that the shape of the orientation distribution strongly affects the rate of the orientation change. At last, the tested closure approximations can lead to good description of some particular cases of orientation distributions. Nonetheless, in general cases, they are not relevant at all so that other strategies would have to be used in order to compute the evolution of fibre orientation.

**Keywords:**
Highly concentrated fibre suspensions
Fibre orientation distribution function
Micro-macro modelling
Computational rheometry

## 1. Introduction

Polymer composites reinforced with short fibre-bundles [1-3] are increasingly used in many industrial processes to produce parts having complex shapes, due to their lightweight and cost-effective properties. The main issue for these materials is to control the fibre content homogeneity and the fibre-bundle orientation induced by the forming processes. Practical applications for structural parts involve highly concentrated suspensions (typical volume fraction of fibres ≈10-30%), where short range bundle-bundle interactions given by relative motion at bundle-bundle contacts are not negligible and can even be regarded as the dominating physical mechanism [4-7].

This observation is the starting point of the approach presented here, which has already been developed and validated in previ-

ous papers [6,7]. In Le Corre et al. [6], a theoretical framework was developed to obtain the mechanical behaviour of a planar network of fibre-bundles using the homogenization technique for periodic discrete structures. In Le Corre et al. [7], efforts aimed at (i) proposing a simple acceptable local bundle-bundle interaction model, (ii) developing numerical simple schemes to mimic real networks of fibre-bundles and to compute their macroscopic behaviour and (iii) comparing the prediction of the micro-macro model with the rheological experiments performed on a non-linear viscous planar suspension, i.e. an industrial Sheet Moulding Compounds (SMC).

Another observation is that numerical simulation works, which aim at predicting the fibre orientation evolution occurring during the flow of this type of industrial concentrated suspensions, are commonly but quite surprisingly based on theories whose domains of validity are *a priori* restricted to dilute or semi-dilute suspensions [8,9]. The models based on this approach are obtained from the well-known Jeffery's equation [10] for a remote fibre in a Newtonian fluid. This model can be modified in order to account for hydrodynamic interactions between fibres by involving the

<footnote>* Corresponding author. Tel.: +33 2 40 37 25 27; fax: +33 2 40 37 25 66.
E-mail address: steven.le-corre@ec-nantes.fr (S. Le Corre).</footnote>
modifications proposed by Tucker and Advani [8]. Usually, this approach leads to describe the evolution of fibre orientation by an equation for the rate of the second-order orientation tensor. This equation has the drawback to include the fourth-order orientation tensor and consequently raises the underlying problem of choosing an appropriate closure approximation to solve the problem. Here the main issue is to investigate the relevance of this type of approach for concentrated suspensions of fibre-bundles.

Hence three ways are followed along this study. Firstly, the predictions of the commonly used models (Jeffery and Folgar Tucker type theories) for fibre orientation evolution are compared with the direct calculation results provided by the presented micro-macro approach. Secondly, the influence of the shape of the fibre orientation distribution on the rate of the orientation evolution is quantified: i.e., a broad spectrum of concentrated suspensions, which may exhibit three different orientation distributions, are numerically tested. Thirdly, as it is known that some closure approximations were developed to work for some characteristic fibre orientation distributions [9], the accuracy of some well-known closure approximations is examined for all considered microstructures.

For that purpose, a non-exhaustive review of commonly used fibre orientation models, based on Jeffery's work, is given in Section 2. In Section 3, a brief presentation of the micro-macro approach, which has been extensively developed in [6,7], is achieved. Some details are given on the numerical technique used to generate the Representative Elementary Volumes (REV's) of the studied suspensions, in particular on the three different considered fibre orientation distributions, i.e., Gaussian-based, uniformly distributed and crenellated. Furthermore, the principle of calculations performed to compute the evolution of the fibre-bundle orientation is explained. Then instantaneous simulation results are presented in Section 4 in order to discuss the influence of the fibre volume fraction and the role of the current orientation state on the rate of the second-order orientation tensor. These results are compared with predictions of the common theories. The same type of comparison is also performed for time-evolution calculations under homogeneous mechanical loading conditions in Section 5.

## 2. State of the art in fibre orientation modelling

### 2.1. Motion of a single fibre immersed in a Newtonian fluid

Jeffery's theory [10] forms the basis of most current approaches to describe fibre suspension rheology and in particular fibre orientation modelling. One of the main contribution of this theory consists in giving an expression for the angular velocity $\omega$ of an ellipsoidal particle in an incompressible infinite Newtonian suspending fluid with no external forces or torques. Later, Hand [11] compared and showed that the Jeffery equation is a special case of a more general equation obtained by Ericksen [12] in his invariant-based theory of anisotropic fluid (the so-called TIF theory) for the evolution of an unit vector $\mathbf{p}$, characterizing a preferred direction at each material point $\mathbf{x}$ of the fluid. If the main axis of revolution of the ellipsoidal particle is described by this unit vector $\mathbf{p}$ and if we consider the cylindrical shape being an approximation to the ellipsoid, the following Jeffery-Ericksen equation may be used to describe the motion of the fibre:

$$ \dot{\mathbf{p}} = \omega \times \mathbf{p} = \mathbf{\Omega} \cdot \mathbf{p} + \lambda (\mathbf{D} \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p})\mathbf{p}), \quad (1) $$

where $\dot{\mathbf{p}}$ denotes the material derivative of $\mathbf{p}$, $\lambda$ is the shape factor of the fibre ($\lambda = (r^2 - 1)/(r^2 + 1)$ with $r = l/d$ the fibre aspect ratio where $l$ and $d$ are respectively the length and the diameter of the fibre), $\mathbf{\Omega}$ the vorticity tensor and $\mathbf{D}$ the strain rate tensor, i.e., the skew-symmetric and symmetric parts of the velocity gradient $\mathbf{L} = \mathbf{gradv}$, respectively, with $\mathbf{v}$ the fluid velocity.

### 2.2. Evolution equation for the orientation distribution function $\psi$

In the case of a suspension of fibres, a reasonable way to describe the orientation of fibres and its evolution consists in adopting a global statistical approach at each material point. For that purpose, the distribution function of orientation $\psi$, as a function of $\mathbf{x}$, $\mathbf{p}$, and $t$, can be used. At each material point $\mathbf{x}$, the elementary probability $dP$ of finding a fibre whose orientation is located in the portion of space delimited by the elementary solid angle $dp$ is written as

$$ dP = \psi(\mathbf{x}, \mathbf{p}, t) dp. \quad (2) $$

The distribution function is $\pi$-periodic: $\psi(\mathbf{x}, \mathbf{p}, t) = \psi(\mathbf{x}, -\mathbf{p}, t)$. It is normalized so that the probability of an orientation on the whole orientation space is unity:

$$ \oint \psi(\mathbf{x}, \mathbf{p}, t) dp = 1. \quad (3) $$

The evolution equation of $\psi$ has the following alternative forms [13]:

$$ \frac{d\psi}{dt} + \psi \mathrm{div}_{\mathbf{p}}(\dot{\mathbf{p}}) = 0, \quad (4) $$

where

$$ \frac{d\psi}{dt} = \frac{\partial \psi}{\partial t} + \nu \cdot \mathbf{grad}_{\mathbf{x}} \psi + \dot{\mathbf{p}} \cdot \mathbf{grad}_{\mathbf{p}} \psi. \quad (5) $$

If the following notation is adopted $D\psi/Dt = (d\psi/dt) - \dot{\mathbf{p}} \cdot \mathbf{grad}_{\mathbf{p}} \psi$, Eq. (4) writes:

$$ \frac{D\psi}{Dt} + \mathrm{div}_{\mathbf{p}}(\psi\dot{\mathbf{p}}) = 0. \quad (6) $$

Following Dupret and Verleye [9], $\dot{\mathbf{p}}$ can be expressed by a constitutive equation, which is a function of $\mathbf{p}$. For instance, the Jeffery equation, whose assumptions are very restrictive, is applicable to describe the evolution of the orientation of fibres contained in dilute suspensions, where hydrodynamic interactions of long and short range between fibres can be considered as negligible.

### 2.3. Evolution equation for the second-order orientation tensor

Determining the orientation state by solving the evolution equation for $\psi$ is quite cumbersome since $\psi(\mathbf{x}, \mathbf{p}, t)$ has six independent variables in the 3D case and this task asks for heavy calculations. This is why orientation tensors were introduced. Orientation tensors are obtained by averaging the $n$ th dyadic products over orientation space, with $\psi$ as the weighting function [14]. Due to the normalization and the periodicity of $\psi$, the odd-order orientation tensors are all zero. Usually only the second and fourth-order orientation tensors are used. They are defined as follows:

$$ \mathbf{A}^{(2)} = \oint \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}) dp, \quad (7) $$

$$ \mathbf{A}^{(4)} = \oint \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}) dp. \quad (8) $$

These tensors only depend on space $\mathbf{x}$ and time $t$ variables. It is then easy to show that the time-evolution equation of the second-order orientation tensor using (4) writes:

$$ \dot{\mathbf{A}}^{(2)} = \oint \{\dot{\mathbf{p}} \otimes \mathbf{p} + \mathbf{p} \otimes \dot{\mathbf{p}}\} \psi(\mathbf{p}) dp. \quad (9) $$

Injecting for instance the form of $\dot{\mathbf{p}}_j$ given by Eq. (1) provides the following equation for $\dot{\mathbf{A}}^{(2)}$:

$$ \dot{\mathbf{A}}^{(2)} = \mathbf{\Omega} \cdot \mathbf{A}^{(2)} - \mathbf{A}^{(2)} \cdot \mathbf{\Omega} + \lambda (\mathbf{A}^{(2)} \cdot \mathbf{D} + \mathbf{D} \cdot \mathbf{A}^{(2)} - 2\mathbf{A}^{(2)} : \mathbf{D}). \quad (10) $$
Notice that, given the imposed macroscopic velocity gradient $L = \mathbf{gradv}$, the above equation implies that $\dot{\mathbf{A}}^{(2)}$ is a linear function of $L$. In particular, we have:

$$ \dot{\mathbf{A}}^{(2)}(-\mathbf{L}) = -\dot{\mathbf{A}}^{(2)}(\mathbf{L}). \quad (11) $$

The evolution equation of the second-order orientation tensor has the drawback to involve the fourth-order orientation tensor $\mathbf{A}^{(4)}$. This difficulty is usually circumvented by defining a closure approximation expressing $\mathbf{A}^{(4)}$ in terms of $\mathbf{A}^{(2)}$. This problem has led to several developments by different authors which first can be found in a paper of Prager [15]. A large variety of closure forms are reviewed in some papers by Advani and Tucker [16] and Dupret and Verleye [9] or Doghri and Tinel [17]. Among the different closure approximations, we will examine in this paper the linear [18,14], the quadratic [19], the hybrid [14] and the natural [9] closure approximations. They are detailed below in their 2D forms, i.e., for fibre orientation contained in a plane $(\mathbf{e}_1, \mathbf{e}_2)$, using indicial notations $(i, j, k, l \in \{1, 2\})$:

* Linear closure approximation [18,14]:

$$ A_{ijkl}^{(4)L} = -\frac{1}{24}(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \frac{1}{6}(A_{ij}\delta_{kl} + A_{ik}\delta_{jl} + A_{il}\delta_{jk} + A_{kl}\delta_{ij} + A_{jl}\delta_{ik} + A_{jk}\delta_{il}), \quad (12) $$

* Quadratic closure approximation [19]:

$$ A_{ijkl}^{(4)Q} = A_{ij}^{(2)}A_{kl}^{(2)}, \quad (13) $$

* Hybrid closure approximation [14]:

$$ A_{ijkl}^{(4)H} = 4 \det \mathbf{A}^{(2)} A_{ijkl}^{(4)linear} + (1 - 4 \det \mathbf{A}^{(2)}) A_{ijkl}^{(4)quadratic}, \quad (14) $$

* Natural closure approximation [9]:

$$ A_{ijkl}^{(4)N} = \frac{1}{3}(A_{ij}A_{kl} + A_{ik}A_{lj} + A_{il}A_{jk}) + \frac{1}{6} \det \mathbf{A}^{(2)}(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{lj} + \delta_{il}\delta_{jk}). \quad (15) $$

## 2.4. Models of hydrodynamic interactions between fibres

In order to take into account fibre-fibre hydrodynamic interactions in semi-dilute suspensions and to provide a constitutive equation for the material derivative $\dot{\mathbf{p}}$, Jeffery's equation can be modified as follows:

$$ \dot{\mathbf{p}} = \dot{\mathbf{p}}_j + \mathbf{q}_d, \quad (16) $$

where $\dot{\mathbf{p}}_j$ is given by (1) and $\mathbf{q}_d$ is a term, which concept in the case of non-Brownian fibre suspension has been first proposed by Folgar and Tucker [13] and later modified and called by Rahnama et al. [20] "orientational dispersion flux". Detailed analyses of the randomization effects of this term on the fibre orientation distribution are given in different studies by Shaqfeh and Koch [21,22] and Rahnama et al. [20] for shear and uniaxial flow kinematics. A general form for this term is proposed by Koch [23] and Petrich et al. [24]:

$$ \mathbf{q}^d = -\mathbf{D}_r \cdot \mathbf{grad}_p \psi, \quad (17) $$

where tensor $\mathbf{D}_r$ is called the orientational diffusion tensor. Different isotropic and anisotropic forms of this diffusion term are proposed. For instance, Folgar and Tucker [13] have first suggested an isotropic diffusion term:

$$ \mathbf{D}_r = \mathbf{D}_r \delta \quad \text{with} \quad \mathbf{D}_r = C_l \dot{\gamma}, \quad (18) $$

where $C_l$ is the interaction coefficient, which needs to be determined from experiments or simulation results (see for instance [25,26] or [27]) and $\dot{\gamma} = \sqrt{2\mathbf{D}:\mathbf{D}}$ the generalized shear strain rate

of the fluid. Some anisotropic diffusion terms have also been proposed. For example, the one given by Fan et al. [28] reads:

$$ \mathbf{D}_r = C_l \dot{\gamma}, \quad (19) $$

where $C_l$ is a second-order symmetric tensor, which components have to be determined experimentally.

In the case of isotropic diffusion, the form of the equation for $\dot{\mathbf{A}}^{(2)}$ is consequently modified as follows [14]:

$$ \dot{\mathbf{A}}^{(2)} = \mathbf{\Omega} \cdot \mathbf{A}^{(2)} - \mathbf{A}^{(2)} \cdot \mathbf{\Omega} + \lambda (\mathbf{A}^{(2)} \cdot \mathbf{D} + \mathbf{D} \cdot \mathbf{A}^{(2)} - 2\mathbf{A}^{(4)} : \mathbf{D}) + 2\mathbf{D}_r(\delta - \alpha\mathbf{A}^{(2)}), \quad (20) $$

with $\alpha = 2$ or 3, respectively in 2D or 3D fibre orientation. To the best of our knowledge, it must be pointed out that all these forms for $\mathbf{D}_r$ break the linearity (11) of $\dot{\mathbf{A}}^{(2)}$ with respect to $L$.

# 3. The present micro-macro approach

This section summarizes the micro-macro model developed in [6,7] to model the rheology of highly concentrated fibre-bundle suspensions with planar bundle orientation.

## 3.1. Micromechanical model and upscaling

The considered suspensions are sheets of thickness $h$ in which slender fibre-bundles of length $l$ ($h < l$) are immersed in an incompressible Newtonian fluid. These bundles form connected networks, i.e., with no isolated bundle or group of bundles. Bundles have elliptical cross-section $\pi d_{max}d_{min}/4$, which main axis $d_{max}(\gg d_{min})$ is contained in the plane of the sheets $\mathcal{P} = (\mathbf{e}_1, \mathbf{e}_2)$. The thickness $h$ being small compared to $l$, it will be assumed that each bundle $b$ has a main and mean orientation along a unit vector $\mathbf{p}_b = \cos\theta_b\mathbf{e}_1 + \sin\theta_b\mathbf{e}_2$ contained in $\mathcal{P}$ [29,30]. This does not mean that bundles are straight: indeed, because of the very high bundle content $f$ ($0.05 < f < 0.3$) and $h < l$, bundles may be initially slightly bent or wavy around their major axis $d_{max}$. Likewise, since $d_{max} \gg d_{min}$, we will (i) neglect bending efforts in the thickness of the sheet, and (ii) neglect the bending of the bundle in $\mathcal{P}$. The last assumption is similar to the one stated in [7,5]. It is supported by experimental evidences. Observations on industrial [3,29] and model [30] concentrated planar bundle suspensions show that the small axis $d_{min}$ of bundles' cross-section in these suspensions are aligned in the $\mathbf{e}_3$-direction. Furthermore, bundles are highly compliant compared with "classical" fibres: they can bend and flatten along $\mathbf{e}_3$ without any significant effort, at least within the investigated strain range. Therefore, this justifies the presented simple micromechanical model for the interactions between bundles where forces in the $\mathbf{e}_3$-direction are neglected. This leads also to the two following consequences:

(a) The motion of a bundle $b$ in $\mathcal{P}$ is simply given by the translational $\tilde{\mathbf{v}}_b$ and angular $\dot{\theta}_b \mathbf{e}_3$ velocities of the centre of mass $G_b$ of the bundle ($\tilde{\mathbf{x}}$ being the projection of $\mathbf{x}$ in $\mathcal{P}$). Moreover, the out-of-plane shearing at the macroscale has not been studied: noting $\mathbf{L} = \mathbf{gradv}$, respectively, with $\mathbf{v}$ the macroscopic velocity gradient of the suspension, only situations for which $L_{i3} = L_{3i} = 0$ ($i = 1, 2$) are explored.

(b) As explained in Le Corre et al. [7], it is assumed that during the deformation of the suspension, contacting bundles interact with forces and local moments induced by the shearing of a small amount of a fluid entrapped in contact zones. As the suspensions under consideration are very concentrated, such deformation mechanisms are supposed to be the main mechanical contributions to the overall stress of the suspension. In this paper, the entrapped fluid is simply considered as Newtonian,
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Modelling of bundle-bundle interactions: (a) view of a connection between two neighbouring fibre-bundles, (b) side view of the sheared zone, (c) top view of the surface of the sheared zone.</figcaption></figure>

its viscosity $\mu$ being that of the suspending fluid. Concerning this last point, very simple assumptions are stated to model the geometry and the mechanics of the entrapped fluid: at a given contact zone $k$ between bundles $b$ and $c$, of centre $K$, such a complex situation is assumed to be equivalent to the shearing of thin prism of height $\alpha$, which in-plane dimensions and orientation depend on $d_{max}$, the orientation of the bundles $\mathbf{p}_b$ and $\mathbf{p}_c$, as well as the relative orientation of the bundles $\Delta\theta_k$ (see Fig. 1). During the relative motion of $b$ and $c$, the entrapped Newtonian fluid is subjected to (i) a simple shear induced by the difference of in-plane translational velocities $\Delta\tilde{v}_k = \tilde{v}_c(K) - \tilde{v}_b(K)$ at point $K$, and (ii) a torsion along $\mathbf{e}_3$ induced by the difference of angular velocities $\Delta\dot{\theta}_k = \dot{\theta}_c - \dot{\theta}_b$. It is then possible to obtain simple estimations of the interaction force $\tilde{\mathbf{f}}_k$ and the local interaction moment $M_k(K)\mathbf{e}_3$ the bundle $c$ exerts on the bundle $b$ [7]:

$$ \tilde{\mathbf{f}}_k = \frac{\mu}{\alpha} \frac{d_{max}^2}{|\sin \Delta\theta_k|} \Delta\tilde{v}_k, \quad M_k(K) = \frac{\pi \mu}{4 \alpha} \frac{d_{max}^3}{|\sin \Delta\theta_k|^2} \Delta\dot{\theta}_k. \quad (21) $$

Therefrom, by neglecting acceleration and external volume forces and moments, by applying the homogenization technique for discrete periodic structures [31–33], and by assuming that the local moments $M_k\mathbf{e}_3$ are small compared to the moments of interaction forces $\tilde{\mathbf{f}}_k$ it can be shown that the self-equilibrium of a Representative Elementary Volume (REV) of the suspensions (of volume $V_{rev}$) can be put in the form [6,7]:

$$ \forall b \in V_{rev}, \quad \left\{ \begin{array}{l} \displaystyle \sum_{k \in C_b} \tilde{\mathbf{f}}_k^{[0]} = \mathbf{0} \\[1.5ex] \displaystyle \sum_{k \in C_b} M_k^{[0]} \mathbf{e}_3 + \sum_{k \in C_b} s_b \mathbf{p}_b \times \tilde{\mathbf{f}}_k^{[0]} = \mathbf{0} \end{array} \right. \quad (22) $$

where $C_b$ is the set of bundles $c$ connected to bundle $b$, and $s_b$ is the local abscissa (from the centre $G_b$ of $b$) along $\mathbf{p}_b$ of the contacting point $K$. In this Eq. (22), $\tilde{\mathbf{f}}_k^{[0]}$ and $M_k^{[0]}$ represent the first-order approximations of $\tilde{\mathbf{f}}_k$ and $M_k\mathbf{e}_3$. They are expressed as follows:

$$ \tilde{\mathbf{f}}_k^{[0]} = \frac{\mu}{\alpha} \frac{d_{max}^2}{|\sin \Delta\theta_k|} \Delta\tilde{v}_k^{[1]}, \quad (23) $$

where

$$ \Delta\tilde{v}_k^{[1]} = \epsilon\tilde{v}_c^{[1]} - \epsilon\tilde{v}_b^{[1]} + \tilde{L} \cdot \tilde{\xi}_k + s_c \dot{\theta}_c^{[0]} \mathbf{e}_3 \times \mathbf{p}_c - s_b \dot{\theta}_b^{[0]} \mathbf{e}_3 \times \mathbf{p}_b, \quad (24) $$

with $\tilde{\xi}_k$ the projection of $G_b G_c$ in $\mathcal{P}$, and

$$ M_k^{[0]} = \frac{\pi \mu}{4 \alpha} \frac{d_{max}^3}{|\sin \Delta\theta_k|^2} (\dot{\theta}_c^{[0]} - \dot{\theta}_b^{[0]}). \quad (25) $$

In these last equations, $\dot{\theta}_b^{[0]}$ is the first-order approximation of $\dot{\theta}_b$, $\epsilon\tilde{v}_b^{[1]}$ is the first-order fluctuation of $\tilde{v}_b$ around its first-order approximation $\tilde{v}_b^{[0]}$, $\tilde{L} = \text{grad}\tilde{v}^{[0]}$ is the macroscopic velocity gradient, which is considered as given and imposed for all fibre-bundles contained in the REV. Please notice that due to the previous kinematical constraints, i.e., incompressibility and $L_{i3} = L_{3i} = 0$ ($i = 1, 2$), this is equivalent to impose $L$. Hence, Eq. (22) represents a system of 3N linear equations with 3N kinematical unknowns, $N$ being the total number of bundles contained in the REV. It is easily shown from (23)–(25) that the unknowns are linear functions of $\tilde{L}$ [6].

Moreover, the homogenization results show the following properties for the suspensions:

(a) the equivalent continuous medium of the studied suspension is a Cauchy medium, whose symmetrical macroscopic stress tensor $\sigma$ is:

$$ \sigma = -p\delta + \tilde{\sigma}^{(f)} \quad \text{with} \quad \tilde{\sigma}^{(f)} = \frac{1}{V_{rev}} \sum_{k=1}^{C} \tilde{\xi}_k \otimes \tilde{\mathbf{f}}_k^{[0]}, \quad (26) $$

where $p$ is the incompressibility pressure, $\delta$ the identity tensor, $V_{rev}$ the volume of the considered REV, $C$ the number of bundle-bundle connections in the REV.

(b) the suspension is a linear viscous anisotropic fluid:

$$ \tilde{\sigma}^{(f)} = \tilde{\sigma}^{(f)}(\tilde{D}, \mu, \text{microstructure}), \quad (27) $$

where $\tilde{D}$ is the strain rate tensor, i.e., the symmetric part of $\tilde{L}$. Such a theoretical result is a consequence of the previous point (a) and of the system of Eqs. (22)–(25).

(c) at the first-order, the rate of change

$$ \dot{\mathbf{p}}_b = \dot{\theta}_b^{[0]} \mathbf{e}_3 \times \mathbf{p}_b = \dot{\theta}_b^{[0]} (\cos \theta_b \mathbf{e}_2 - \sin \theta_b \mathbf{e}_1) \quad (28) $$
of any orientation vector $\mathbf{p}_b$ is a linear function of $\tilde{L}$:

$$ \left. \begin{array}{l} \forall i, j \in \{1, 2\} \\ \forall L_{ij} \in \mathbb{R} \\ \forall b \in V_{\text{rev}} \end{array} \right\} \dot{\mathbf{p}}_b (L_{ij} \mathbf{e}_i \otimes \mathbf{e}_j) = L_{ij} \dot{\mathbf{p}}_b (\mathbf{e}_i \otimes \mathbf{e}_j), \qquad (29) $$

where the Einstein convention has been used for the summation on indices i and j. This also a direct consequence of the system of Eqs. (22)–(25).

(d) Consequently, the rate of the corresponding distribution function $\psi$ and of the n th-order orientation tensors $\mathbf{A}^{(n)}$ are also linear functions of $\mathbf{L}$. In particular, property (11) is also verified.

## 3.2. Microstructure generation

In order to mimic real fibrous microstructures of industrial polymer composites such as Sheet Molding Compounds (SMC) or Glass Mat Thermoplastics (GMT), a very simple deterministic technique directly inspired from the statistical tube model [34–36] has been developed [7]. Following this tube model, the number of bundle-bundle contacts per unit volume should be a quadratic function of the volume fraction $f$ of bundles. Including the tube model in the previous homogenization approach [7] induces that the macroscopic extra stress tensor, (e.g., the term $\tilde{\sigma}^{(f)}$ of the Eq. (26)) is proportional to the number of contacts in the suspension, and thus evolves as a quadratic function of the fibre volume fraction $f$. Experimental results [3,30] depict also such quadratic evolution of stress levels with $f$: this reinforces the relevance of the tube model despite its simplicity.

Therefore, this generation process will be used in the following sections to study the evolution of the orientation of the bundles during mechanical loading. Briefly, $N$ bundles (elementary volumes $\pi d_{max}d_{min}l/4$) which orientation vectors $\mathbf{p}_b$ are contained in $\mathcal{P}$ are generated inside a 3D REV of volume $V_{\text{rev}} = d \times d \times h$ (see Fig. 2, $h$ being the thickness of the REV). The centre of mass $G_b$ of each bundle has a random position $(x_b, y_b, z_b)$ inside the volume of the REV.

The in-plane dimension $d$ of the REV must be chosen larger than the bundle length $l$, in order to avoid the possibility to generate a continuous fibre-bundle. Practically, $d$ was set to a length 1% higher than $l$, in order to enhance the computation time. Within the bundle content range for which the micro-macro approach was initially developed (>8–10% for the considered bundle aspect ratio), i.e., when the number of bundle-bundle contacts per bundle is greater or equal to 2 (largely above the percolation threshold), a preliminary study has proved that numerical results were not affected by this choice of $d$. At lower fibre-bundle content, both the physics at the bundle scale and the adopted generation procedure may become questionable so that the size of the REV should be reconsidered and increased (see for instance Section 4.1).

The generation process permits to obtain three different shapes for the fibre orientation distribution of bundle angles $\theta_b$ of the studied planar networks, i.e., uniform, Gaussian-based and crenellated, the two latter ones being centred on $\theta = 0$. They are shown in Fig. 2(a-c). The characteristics of the uniform distribution are that its amplitude is chosen and that there are no privileged classes of fibre angles $\theta_b$ (cf. Fig. 2(a)). The Gaussian-based distribution consists first in starting from a Gaussian distribution of orientation angles $\theta_b$ for fibres, and then in using the $\pi$-periodicity of $\psi$ with respect to $\theta_b$ to normalize the resulting histograms (cf. Fig. 2(b)). The crenellated distribution is detailed hereafter. It exhibits a peak, which intensity and amplitude can be varied. More precisely, given the number $N$ of fibres in the REV, two orientation amplitudes $\alpha_1$ and $\alpha_2$ are imposed to two populations of $N_1$ and $N_2 = N - N_1$ fibres, respectively. A random distribution is set for both popula-

tions of fibres, leading to superimpose two random distributions between $-\alpha_1/2$ and $+\alpha_1/2$ and between $-\alpha_2/2$ and $+\alpha_2/2$, respectively (fibrous microstructures are oriented in the $\mathbf{e}_1$-direction). As depicted in the simplified graph of Fig. 2(c), fibre orientation in generated REV's can be characterized by $\alpha_1$, the base orientation amplitude, $r_\alpha = \alpha_2/\alpha_1$, the amplitude ratio and $r_N = N_2/N$, the fraction of fibres in the peak exhibited by this type of orientation distribution. Notice that at fixed values of $N$ and $\alpha_1$, the increase of $r_N$ is accompanied with an increase of $r_\alpha$.

In order to detect bundle-bundle contacts, a control volume $V_b$ is assigned to each bundle $b$: every bundle $c$ which centre-line intersects $V_b$ is added to the connectivity set of bundle $b$. The control volume $V_b$ is chosen as a rectangular box with dimensions $l \times d_{max} \times \delta^* d_{min}$, where $\delta^*$ is a dimensionless parameter [7]. In the following, the generated bundle networks will be characterized by their bundle volume fraction $f$, and by the second $\tilde{\mathbf{A}}^{(2)}$ and fourth $\tilde{\mathbf{A}}^{(4)}$ 2D-orientation tensors. These macroscopic descriptors of the microstructure are respectively defined in a discrete way by:

$$ \begin{gathered} f = N \frac{\pi d_{max} d_{min} l}{4 V_{\text{rev}}}, \\ \tilde{\mathbf{A}}^{(2)} = \frac{1}{N} \sum_{b=1}^{N} \tilde{\mathbf{p}}_b \otimes \tilde{\mathbf{p}}_b, \quad \tilde{\mathbf{A}}^{(4)} = \frac{1}{N} \sum_{b=1}^{N} \tilde{\mathbf{p}}_b \otimes \tilde{\mathbf{p}}_b \otimes \tilde{\mathbf{p}}_b \otimes \tilde{\mathbf{p}}_b. \end{gathered} \qquad (30) $$

## 3.3. Computational rheometry

Two types of calculations will be presented in the following:

(1) **Instantaneous mechanical responses of the generated REV's**

—Thanks to the linearity property (29), the mechanical response of a microstructure at a given time can be entirely determined imposing on the studied REV's four independent macroscopic velocity gradients, i.e., $\mathbf{e}_i \otimes \mathbf{e}_j(i, j \in \{1, 2\})$. In practice, given one of these four macroscopic velocity gradients, this consists in inverting the linear system (22) in order to compute (i) the unknown kinematical fields $\epsilon \tilde{\mathbf{v}}_b^{[1]}$ and $\dot{\tilde{\mathbf{v}}}_b^{[0]}$, (ii) the overall stress tensor of the suspension $\sigma$ (according to (26)) and (iii) the rate $\dot{\tilde{\mathbf{A}}}^{(2)}$ of the orientation tensor $\tilde{\mathbf{A}}^{(2)}$ defined as (according to (28) and (30)):

$$ \dot{\tilde{\mathbf{A}}}^{(2)} = \frac{1}{N} \sum_{b=1}^{N} (\dot{\tilde{\mathbf{p}}}_b \otimes \tilde{\mathbf{p}}_b + \tilde{\mathbf{p}}_b \otimes \dot{\tilde{\mathbf{p}}}_b). \qquad (31) $$

(2) **Time-evolution of microstructures**

—In order to determine the time-evolution of a REV during a given mechanical loading, both a Lagrangian point of view and an explicit time-integration scheme have been adopted. More precisely, after each “snapshot calculation” performed at a given time $t$, the position $\mathbf{x}_b^{t+\delta t}$ of the centre of mass $G_b$ as well as the orientation $\mathbf{p}_b^{t+\delta t}$ of a bundle $b$ at time $t + \delta t$ are updated from the knowledge of the position $\mathbf{x}_b^t$, the orientation $\mathbf{p}_b^t$ and the computed velocity fields $\epsilon \mathbf{v}_b^{[1]}$ and $\dot{\mathbf{v}}_b^{[0]} \mathbf{e}_3$ at $t$:

$$ \left\{ \begin{aligned} \mathbf{x}_b^{t+\delta t} &= \mathbf{x}_b^t + (\mathbf{L} \cdot \mathbf{x}_b^t + \epsilon \mathbf{v}_b^{[1]}) \delta t \\ \mathbf{p}_b^{t+\delta t} &= \mathbf{p}_b^t + (\dot{\mathbf{v}}_b^{[0]} \mathbf{e}_3 \times \mathbf{p}_b^t) \delta t. \end{aligned} \right. \qquad (32) $$

Therewith, the connectivity of the updated REV is achieved following the method given in Section 3.2, and a new time iteration can be achieved. A preliminary study has been performed in order to determine the time-step values $\delta t$ that give reproducible and time-step insensitive simulations.
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** (a) Uniform fibre orientation distribution ranging between -60° and 60°, (b) Gaussian-based fibre orientation distribution, (c) crenellated fibre orientation distribution, with their corresponding second-order orientation tensors $\tilde{A}^{(2)}$, (d) example of a generated REV ($V^{(s)} = 25 \text{ mm} \times 25 \text{ mm} \times 2.5 \text{ mm}$) containing 437 bundles of length 25 mm ($f = 0.188$) (symbols ○ represent the connections location, lines represent the mean center-line of the bundles). (e) Same REV after flowing under a unit plane strain deformation until an axial Hencky strain $\epsilon_{33} = \ln(h/h_0) = -1$.</figcaption></figure>

In order (i) to underline the role of both the fibrous microstructure (bundle content and orientation) and the loading path on the macroscopic response of the suspension, and (ii) to test the capability of usual fibre orientation models in the concentrated regime, numerical simulations presented below were performed with the following microstructural parameters which values are typical of those encountered in polymer composites (SMC or GMT):

* the dimensions of the bundles were $l = 25 \text{ mm}$, $d_{max} = 0.6 \text{ mm}$, $d_{min} = 0.06 \text{ mm}$, and the REV's are $25 \text{ mm} \times 25 \text{ mm} \times 2.5 \text{ mm}$ rectangular boxes [7],

* bundle networks were generated with the procedure described in Section 3.2($\tilde{A}_{12}^{(2)} \approx 0$). The principal orientation of the networks $\tilde{A}_{11}^{(2)}$ ranged between 0 and 1, the volume fraction of bundles $f$ between 0.05 and 0.3,

* the thickness $\alpha$ of the equivalent sheared contact prisms was fixed to $2 \mu\text{m}$, the dimensionless parameter $\delta^*$ to 3 [7].

The generated REV's were subjected to two mechanical loading conditions:

* a unit plane strain compression (l) in ($e_1, e_3$), i.e., a 1D elongational flow in $\mathcal{P}$:
$$L_i = e_1 \otimes e_1 - e_3 \otimes e_3, \qquad (33)$$

* a unit pure shear (II) in ($e_1, e_2$), i.e., an irrotational shear flow in $\mathcal{P}$:

$$L_{II} = \frac{1}{2}(e_1 \otimes e_2 + e_2 \otimes e_1). \qquad (34)$$

As an example, Fig. 2(b) gives the deformed shape of the REV plotted in Fig. 2(a) after a unit plane strain compression up to an axial Hencky strain $\epsilon_{33} = \ln(h/h_0) = -1$.

## 4. Instantaneous responses

The first presented simulation results are snapshots in time of the instantaneous evolution of fibre-bundle networks. The following sections aim at describing the influence of the fibre volume fraction or the influence of current fibre orientation on the rate of orientation $\dot{\mathbf{p}}_b$ and on the orientation rate tensor $\dot{\mathbf{A}}^{(2)}$. This is first performed using micro-macro modelling predictions (cf. Sections 4.1 and 4.2). Then a comparison with usual approaches (cf. Section 2) is achieved in Sections 4.3–4.5.

### 4.1. Influence of the fibre volume fraction on the orientation rate

Fig. 3 (a) gives the evolution of the components $\dot{\mathbf{A}}_{11}^{(2)}$ and $\dot{\mathbf{A}}_{12}^{(2)}$ of the tensor $\dot{\mathbf{A}}^{(2)}$ calculated for nearly isotropic Gaussian-based

microstructures ($\dot{\mathbf{A}}_{11}^{(2)} \approx 0.5$) and subjected to a 1D elongational flow. Fig. 3(b) gives the evolutions of the same components for highly oriented microstructures ($\dot{\mathbf{A}}_{11}^{(2)} \approx 0.8$) tested for the same flow kinematics. Whatever the studied microstructure (isotropic or oriented), we can notice that the $\dot{\mathbf{A}}_{12}^{(2)}$ components are very low compared to the $\dot{\mathbf{A}}_{11}^{(2)}$ components.

* For $f > 0.05$, it appears that the fibre volume fraction has not an important influence on the results.

* For $f \le 0.05$, the same observation can be made despite a larger scattering of the results. The origin of this scattering is mainly due to the low number of bundles that are present in REV's, moreover some remote bundles or remote clusters of bundles can be found in the generated REV's. These elements are not accounted for the calculations of the macroscopic properties of the generated REV's [7]. Typical sizes of these clusters can be of the same order as the sizes of the REV's so that the size of the REV's may be inappropriate for such situations. To obtain reliable results two alternatives can be followed: the first one consists in increasing the size of the generated REV's, the second one consists in performing a large number of calculations on REV's of limited sizes in order to get representative average values of the studied properties [37,38]. Here, this last option was chosen. It has also to be pointed out that remote bundles or clusters could have an impact on the rate of the orientation tensor in this range of rather low fibre volume fractions. In this particular case, the micro-macro approach

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Influence of the fibre volume fraction $f$ on the rate of the components of the orientation tensor $\dot{\mathbf{A}}^{(2)}$ for 1D elongational flows (a and b) and for pure shear flows (c and d) in $\mathcal{P}$- (a) and (c): nearly isotropic REV's ($\dot{\mathbf{A}}_{11}^{(2)} \approx 0.5$) - (b) and (d): highly oriented REV's along $e_1$ ($\dot{\mathbf{A}}_{11}^{(2)} \approx 0.8$). Calculations are performed using Gaussian-based fibre orientation distributions.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Components $\dot{\tilde{p}}_i$ of the orientation evolution vector $\dot{\tilde{p}}_b$ for each fibre-bundle of REV's having a number of fibre-bundles $N = 600$. Comparison between micro-macro (symbols) and "Jeffery"(lines) predictions for 1D elongational flows in $\mathcal{P}$ along $\mathbf{e}_1$. Calculations are performed using uniform fibre orientation distributions. (a) Nearly isotropic - $f = 10\%$, (b) highly oriented - $f = 10\%$, (c) nearly isotropic - $f = 25\%$, (d) highly oriented - $f = 25\%$.</figcaption></figure>

should be improved by accounting for bundle-fluid interactions in order to incorporate the influence of these remote elements or domains in the description of the rheological behaviour of the suspension.

* When these microstructures are subjected to a pure shear kinematics (Fig. 3(c) and (d)), the same trends are observed.

* Fig. 4 (a–d) shows the components of $\tilde{\mathbf{p}}_b$ obtained by the micro-macro approach for each bundle of nearly isotropic and very oriented cells subjected to 1D elongational flows. Two fibre-bundle fractions ($f = 10\%$ and $25\%$) are considered. These figures also underline the strong influence of the global orientation whatever the fibre fraction $f$. For a given orientation, clouds of points represented in these figures are close when comparing the two tested fibre-bundle fractions. Nonetheless, it has to be noticed that the dispersion of results decreases when increasing $f$: i.e., the clouds of points is more compact for high bundle fractions. These observations have to be related to the decrease of the dispersion of the components of $\tilde{\mathbf{A}}^{(2)}$ of Fig. 3(a–d).

* Thus, given the very weak influence of the fibre volume fraction $f$ on the components of $\tilde{\mathbf{A}}^{(2)}$ for elongational and shear flows for

$f > 0.05$, we will focus the rest of our study on differently oriented microstructures with an identical fibre content of 0.2.

## 4.2. Influence of the current orientation of fibres on the orientation rate

First results of Fig. 5 (a) and (b) are obtained from calculations performed using cells where the fibre orientation distribution is Gaussian-based. In Fig. 5(a), the components of $\tilde{\mathbf{A}}^{(2)}$ are given as a function of the current orientation $\tilde{\mathbf{A}}_{11}^{(2)}$ for cells having different fibre orientations subjected to an elongational kinematics ($f = 0.2$). Micro-macro predictions show that the orientation rate in the flow direction (given by $\tilde{\mathbf{A}}_{11}^{(2)}$) strongly depends on the current fibre orientation, which is given here by $\tilde{\mathbf{A}}_{11}^{(2)}$. The maximum of fibre orientation evolution in the direction of the flow is obtained for nearly isotropic structures ($\tilde{\mathbf{A}}_{11}^{(2)} \approx 0.5$) and this decreases with the fibre orientation intensity, either in the $\mathbf{e}_1$-direction or in the $\mathbf{e}_2$-direction. It is interesting to notice that the evolution of $\tilde{\mathbf{A}}_{11}^{(2)}$ is not exactly symmetrical with respect to $\tilde{\mathbf{A}}_{11}^{(2)} = 0.5$: the decrease of
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Influence of the orientation intensity $\tilde{\mathbf{A}}_{11}^{(2)}$ on the components of the orientation tensor rate $\dot{\tilde{\mathbf{A}}}^{(2)}$: (a) for 1D elongational flows in $\mathbf{P}$ along $\mathbf{e}_1$, (b) for pure shear in the case of Gaussian-based fibre orientation distributions, (c) and (d) are respectively the same as (a) and (b) in the case of uniform fibre orientation distributions. Comparison between micro-macro and “Jeffery” predictions.</figcaption></figure>

$\dot{\tilde{\mathbf{A}}}_{11}^{(2)}$ is higher when $\tilde{\mathbf{A}}_{11}^{(2)} < 0.5$ than when $\tilde{\mathbf{A}}_{11}^{(2)} \ge 0.5$. Notice that $\tilde{\mathbf{A}}_{11}^{(2)}$ has nearly no influence on $\dot{\tilde{\mathbf{A}}}_{12}^{(2)}$. Even if it is never exactly zero, it remains small compared to $\dot{\tilde{\mathbf{A}}}_{11}^{(2)}$. This can simply be attributed to the random generation process where the average alignment of fibres along $\mathbf{e}_1$ is never perfect.

Fig. 5(b) represents the same kind of results for oriented REV's subjected to pure shear kinematics. In this case, the first principal value $\dot{\tilde{\mathbf{A}}}_1^{(2)}$ and the value of the angle between the major principal vector $\mathbf{e}_1$ and the direction $\mathbf{e}_1$ of $\dot{\tilde{\mathbf{A}}}^{(2)}$ are represented with respect to $\tilde{\mathbf{A}}_{11}^{(2)}$. The micro-macro predictions show that this angle is close to $\pi/4$ whatever the current fibre orientation given by $\tilde{\mathbf{A}}_{11}^{(2)}$. Hence, even if the subjected kinematical loading is irrotational, the fibrous microstructure rotates. The maximum of $\dot{\tilde{\mathbf{A}}}_1^{(2)}$ is obtained for cells highly oriented either in the $\mathbf{e}_1$-direction or in the $\mathbf{e}_2$-direction. Moreover, this evolution appears to be symmetrical with respect to $\tilde{\mathbf{A}}_{11}^{(2)} = 0.5$ in this case.

Fig. 5(c) and (d) give respectively the same type of results as Fig. 5(a) and (b) for calculations performed using uniform fibre orientation distributions. General trends are quite close to the previous figures. Nevertheless, the evolution of $\dot{\tilde{\mathbf{A}}}_{11}^{(2)}$ (respectively $\tilde{\mathbf{A}}_1^{(2)}$) exhibits a more complex “wavy” shape with two local maxima

(respectively minima) for the elongational kinematics (respectively the pure shear) than in the case of the Gaussian-based calculations. Moreover, the obtained rates can differ quite largely too, e.g., $\dot{\tilde{\mathbf{A}}}_{11}^{(2)} \approx 0.3\,\text{s}^{-1}$ at $\tilde{\mathbf{A}}_{11}^{(2)} = 0.35$ for the uniform distributions whereas $\dot{\tilde{\mathbf{A}}}_{11}^{(2)} \approx 0.24\,\text{s}^{-1}$ for the Gaussian-based case and 1D elongational flows. Thus, it appears that the shape of the fibre orientation distribution is preponderant on the rate of the evolution of the orientation of fibres whatever the tested kinematics.

This important observation is confirmed when using crenellated fibre orientation distributions as depicted in Fig. 6 (a) for elongational kinematics. This figure indicates that the rate of the components of the second-order orientation tensor will “slow down” as the fraction fibres $r_N$ of fibres will decrease, i.e., when the orientation distributions exhibit sharp orientation peaks. As evident from the graph of Fig. 6(c), for cells having a similar average value $\langle \tilde{\mathbf{A}}_{11}^{(2)} \rangle$ of the orientation tensor $(\tilde{\mathbf{A}}_{11}^{(2)} \in [0.64, 0.68])$, a very large scattering of the orientation rates can be induced by changing the orientation parameters $r_N$ and $r_\alpha$: $0.17\,\text{s}^{-1} \langle \tilde{\mathbf{A}}_{11}^{(2)} \rangle < 0.33\,\text{s}^{-1}$. This figure shows that the lower the number fraction $r_N$, i.e., the sharper the peak shape, then the lower the orientation rate $\tilde{\mathbf{A}}_{11}^{(2)}$. This trend is systematically observed for other orientation states $\langle \tilde{\mathbf{A}}_{11}^{(2)} \rangle > 0.5$. It can also be noticed that calculations carried out with crenel-
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Influence of the orientation intensity $\tilde{\mathbf{A}}_{11}^{(2)}$ on $\tilde{\mathbf{A}}_{11}^{(2)}$, for 1D elongational flows in $\mathcal{P}$ along $\mathbf{e}_1$. (a) Comparison of the micro-macro predictions in the cases of uniform, Gaussian-based or crenellated fibre orientation distributions with the predictions of the hybrid and natural closure approximations, (b) same type of comparison when using this time the Jeffery predictions instead of the micro-macro ones. Remark: only polynomial best fits of the simulation results are given to make the figures readable. (c) Evolution of $\tilde{\mathbf{A}}_{11}^{(2)}$ for the micro-macro approach with respect to $r_N$ ($\alpha_1 = \pi, 0 < r_\alpha < 1$, REV's with a similar average orientation such that $\tilde{\mathbf{A}}_{11}^{(2)} \in [0.64, 0.68]$).</figcaption></figure>

lated fibre orientation distributions with $r_N > 0.8$ are naturally close to the results obtained for uniform orientation distributions. A simple rule therefore can be established: fibre orientation distributions which exhibit some peaks, i.e., like the Gaussian-based or crenellated orientation distributions, will give slower rates for the orientation evolution than the uniform distribution. The relevance of a particular type of fibre orientation distributions is investigated in Section 5.

### 4.3. Influence of the current orientation of fibres on the orientation rate tensor using “Jeffery” predictions

Results given in Fig. 5 were also obtained for simulations called “Jeffery” simulations. The “Jeffery” simulations consisted in calculating the rate of orientation change of each generated single fibre of the tested cells adopting Jeffery's Eq. (1) and then using these data to calculate the components of the orientation rate tensor $\dot{\tilde{\mathbf{A}}}^{(2)}$ according to Eq. (31). Notice that another equivalent approach would have consisted in performing a calculation of Eq. (10), which is only the rewritten Jeffery equation for the rate of the second-order orientation tensor, using a discrete evaluation of $\tilde{\mathbf{A}}^{(2)}$ and $\tilde{\mathbf{A}}^{(4)}$. It was checked that both ways of calculation gave similar results.

In a sense, our approach could be seen as using an exact closure approximation for $\tilde{\mathbf{A}}^{(4)}$ in the expression for $\dot{\tilde{\mathbf{A}}}^{(2)}$ given by Eq. (10).

When comparing this approach to the micro-macro computations one can check that the predicted trends are, at a first sight, very similar, which might seem surprising. This observation holds whatever the tested orientation Gaussian-based and uniform distributions. For the pure shear case, the shape of the both predictions are very close all over the orientation range. For the elongational flow, it is interesting to notice that both approaches do not match exactly. If the Jeffery prediction gives a symmetrical evolution with respect to $\tilde{\mathbf{A}}_{11}^{(2)} = 0.5$, this is no more the case of results predicted by the micro-macro approach. This is particularly evident for rather oriented networks ($\tilde{\mathbf{A}}_{11}^{(2)} > 0.7$) along the flow direction. This difference should be further studied, but could possibly be attributed to the pronounced anisotropy of the fibrous medium in these situations: e.g., the Jeffery equation considers the case of the orientation evolution of a remote fibre immersed in a isotropic viscous medium, but here the situation is largely different and one would rather have to consider the case of a fibre immersed in an anisotropic medium [22]. It should be also pointed out that the influence of local interaction moments might increase, as shown in Le Corre et al. [7] for these particularly oriented microstructures. They might possibly
become no more negligible at the macroscopic level. Accounting for their effects in the micro-macro approach could be possible: this would lead to an equivalent description of the rheological behaviour of the suspension, which would be a general Cosserat medium, as shown theoretically in [6].

Fig. 4(a–d) allows comparing the components of $\dot{\mathbf{p}}_b$ given for individual cells by the micro-macro and the Jeffery approaches. Predicted shapes of the clouds of points given by the micro-macro approach are quite close to the Jeffery prediction, whatever the fibre-bundle fractions and orientations, but the micro-macro one depicts a more complex situation than the Jeffery one. Indeed, all results of the micro-macro approach appear to be dispersed around the lines given by using the Jeffery theory. This is especially the case of nearly isotropic microstructures, but it seems that bigger differences can be observed when the microstructures are highly oriented. This latter result has to be linked with the previous observation on the evolution of the components of $\dot{\mathbf{A}}^{(2)}$.

Jeffery predictions are also close to the micro-macro ones for crenellated orientation distributions as it appears when comparing Fig. 6(a) and (b). Thus, the relative accuracy of the Jeffery prediction with respect to the micro-macro one (except for highly oriented microstructures) seems to be a general rule.

### 4.4. Influence of the current orientation of fibres on the orientation rate tensor using closure approximations

Results of Fig. 7 are obtained for Gaussian-based fibre orientation distributions. They enable to compare the predictions for the evolutions of the components of the orientation rate tensor $\dot{\mathbf{A}}^{(2)}$ given by the micro-macro approach with the predictions of usual modelling approaches using closure approximations in Eq. (10) for 1D elongational and pure shear kinematics. As in Section 4.3, the Jeffery approach is also placed in those figures for comparison. Whatever the tested kinematics, the following points are worth being noted:

* For isotropic states of fibre orientation ($\tilde{\mathbf{A}}_{11}^{(2)}$ close to 0.5), the linear closure approximation gives quite good predictions when compared to the micro-macro ones. This result was expected as the linear closure approximation is exact for isotropic orientation states if every fibre rotates according to Jeffery's equation (see

e.g. [9]). Nonetheless, the predicted orientation rate appears to be insensitive to current orientation intensity, which is a serious drawback. For instance, this approximation does not predict that the orientation decreases for very oriented states and therefore seems to seriously fails.

* At the opposite, the quadratic closure approximation largely overestimates the $\tilde{\mathbf{A}}_{12}^{(2)}$ and $\tilde{\mathbf{A}}_{11}^{(2)}$ for both tested kinematics and isotropic orientation and is therefore adapted to very oriented states only. As previously, this latter result was also expected as the quadratic closure is exact for any aligned orientation states if every fibre rotates according to Jeffery's equation (see e.g. [9]).

* The hybrid and natural closure approximations follow correctly the trends given by the micro-macro approach and by Jeffery's type calculations. Their predictions are good for isotropic and very oriented states and a certain discrepancy is observed for intermediate regimes, depending on the mechanical loading path.

* Whatever the tested kinematics, usual modelling approaches predict that the evolution of the components of $\dot{\mathbf{A}}^{(2)}$ is symmetrical with respect to the isotropic fibre orientation state. If this symmetry is observed for the micro-macro approach in the case of the pure shear kinematics, this is not the case, as we have previously noticed it, for the 1D elongational kinematics.

Likewise, predictions for the components of $\dot{\mathbf{A}}^{(2)}$ of the hybrid and natural closure approximations are compared to the micro-macro and Jeffery results for crenellated and uniform fibre orientation distributions in Fig. 6 for the 1D elongational kinematics. Whatever the tested orientation distributions, none of those closure approximation gives results corresponding to the micro-macro or Jeffery results over the whole range investigated for $\tilde{\mathbf{A}}^{(2)}$. Nevertheless, it can be observed that these two closure approximations give trends that belong to the large spectrum of possible evolution rates that can be reached for all the possible orientation distributions and orientation intensities.

To summarize, no closure approximation compares with the micro-macro approach like the Jeffery-type approach, whatever the tested orientation distributions and the kinematics, which are not exhaustively presented here. It should be noticed that the computation of the Jeffery approach induces an exact calculation of $\tilde{\mathbf{A}}^{(4)}$ in a discrete way (see also the remark given at the beginning of the previous section). Thus, the Jeffery approach, which

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Influence of the orientation intensity $\tilde{\mathbf{A}}_{11}^{(2)}$ on the components of the orientation tensor rate $\dot{\mathbf{A}}^{(2)}$: (a) for 1D elongational flows in $\mathcal{P}$ along $e_1$, (b) for pure shear. Calculations are performed using Gaussian-based fibre orientation distributions. Comparison between micro-macro predictions and closure approximations predictions.</figcaption></figure>
adopts a mean field assumption for the macroscopic velocity field, seems to be the most relevant when compared to the micro-macro one. Nevertheless, none of those approaches explains the observed discrepancy for elongational flows of oriented fibre networks ($\tilde{A}_{11}^{(2)} > 0.7$).

### 4.5. Hydrodynamic diffusion terms

As mentioned in Section 2.4, corrective terms, known as hydrodynamic diffusion terms, have been added by several authors in order to account for the fibre-fibre interactions, in the case semidilute of suspensions [25,27]. In this section, we compare the predictions of the micro-macro approach with those of the Jeffery approach using a similar technique to compute $\tilde{A}^{(2)}$ as in Sections 4.3 and 4.4, but adding one term, which accounts for the hydrodynamic diffusion. All tested microstructures of this section have Gaussian-based orientation distributions. One of the simplest diffusion model proposed in the literature by Phan-Thien et al.'s (PT) is chosen here. This diffusion term implies to use an interaction coefficient, which is written as [27]:

$$C_I = 0.03(1 - \exp(-0.224fr)). \quad (35)$$

Here, this coefficient is equal to $C_I = 0.0251$ for $f = 0.2$. It is important to make the following remarks concerning this model.

* On the contrary to this model, it has to be noticed that no influence of the fibre volume fraction $f$ ($f \ge 0.05$) is revealed by the micro-macro approach.

* Adding this term to the equation for the evolution of $\tilde{A}^{(2)}$ leads to break its linearity (as previously mentioned in Section 2) on the contrary to the simplest Jeffery or micro-macro predictions.

As clear from Fig. 8, the interaction coefficient influence is not consistent with the trends obtained by the micro-macro approach for the studied concentrated suspensions. In Fig. 8(a), the orientation rate for $\tilde{A}_{11}^{(2)}$ is overestimated for low orientation whereas it is underestimated for high orientation for 1D elongational flow cases. Thus, it seems that the increase of the orientation rate observed for highly oriented microstructures in the direction of the 1D elongation cannot be described by this type of interaction coefficient. For pure shear loading (cf. Fig. 8(b)) the evolution of $\tilde{A}_{12}^{(2)}$ is well

predicted, but the evolutions of $\tilde{A}_{11}^{(2)}$, which are obtained by the micro-macro or PT approaches are drastically different. In the case of the PT model, $\tilde{A}_{11}^{(2)}$ decreases as $\tilde{A}_{11}^{(2)}$ increases, whereas it remains close to zero in the case of the micro-macro approach. This observation means that both predictions will give completely different principal orientation axes for the 2D second-order orientation tensor $\tilde{A}^{(2)}$.

### 5. Time-evolution of fibre orientation

The evolution of fibrous microstructures subjected for 1 s to a 1D elongational flow (constant plane strain compression at $1 \text{ s}^{-1}$, see (33)) is now analyzed. Notice that simulations were performed with a reliable time-step of $10^{-2} \text{ s}$. Also notice that such flow conditions correspond to a REV height reduction of about 65 %, which is close to typical processing conditions. Two different methods were used to compute the evolution of fibre orientations and spatial positions within the tested REV's:

* The first method consisted in performing “direct” calculations, by integrating step by step with Eq. (32) the fibre orientation and spatial position of each fibre contained in REV's. This was done by using the micro-macro and the Jeffery-type approaches. For the first approach, the translation $\epsilon v_b^{[1]}$ and rotation $\dot{\theta}_b^{[0]}$ velocity fields were obtained by solving the localization problem (22). For the second approach, the rotation velocities $\dot{\theta}_b^{[0]}$ were computed with Jeffery's Eq. (1), and the $\epsilon v_b^{[1]}$'s were set to 0 (affine motion of fibres). The initial fibrous microstructures displayed random fibre positions and three Gaussian-based distributions of fibre orientation, aligned in the $e_1$-direction and which initial orientation tensors were such that $\tilde{A}_{11}^{(2)} = 0.25, 0.5$ and $0.75$, respectively. All results, that were averaged over 20 REV's with the same microstructural parameters, are presented in Fig. 9, which gives the evolutions of $\tilde{A}_{11}^{(2)}, \tilde{A}_{11}^{(2)}, \tilde{A}_{12}^{(2)}$ and $\tilde{A}_{12}^{(2)}$ with respect to time (a) and $\tilde{A}_{11}^{(2)}$ (b), respectively.

* The second method, also carried out with the same initial microstructures, assumed that the positions of fibres always remained randomly spatially distributed in the REV's, and that the orientations of fibres always exhibited Gaussian-based distributions. In order to compute the evolution of fibre orientation, the values of $\tilde{A}_{11}^{(2)}$ were obtained from polynomial fits of the data

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Influence of the orientation intensity $\tilde{A}_{11}^{(2)}$ on the components of the orientation tensor rate $\dot{\tilde{A}}^{(2)}$: (a) for 1D elongational flows along $e_1$, (b) for pure shear. Calculations are performed using Gaussian-based fibre orientation distributions. Comparison with the predictions of Phan-Thien et al.'s interaction model [27].</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Average orientation evolution in time for 20 different cells with initially Gaussian-based distributions: planar compression test along *e*₁ with *D*₁₁ = −*D*₃₃ = 1 s⁻¹, *N* = 500 fibres: (*a*ᵢ) Ā⁽²⁾ components versus time, (*b*ᵢ) Ā⁽²⁾ components versus Ā₁₁⁽²⁾ - initial value of Ā₁₁⁽²⁾: 0.25 (x₁), 0.50 (x₂), 0.75 (x₃).</figcaption></figure>
plotted in Figs. 7(a) and 8(a). Results deduced from this method have also been reported in Fig. 9 (the so-called “fit” results).

* Likewise, Fig. 11 displays results obtained by the first method for the micro-macro and Jeffery calculations. Results called “hybrid” and “natural” are given by adopting the second method described in the previous point. The values of $\tilde{A}_{11}^{(2)}$ are fitted from the data of Fig. 7 for the natural and hybrid closure approximations. Here, an average initial orientation of the tested microstructures is taken into account, i.e., $\tilde{A}_{11}^{(2)} = 0.5$.

* At last, Fig. 12 shows similar results for initially crenellated orientation distributions with approximately the same value of $\tilde{A}_{11}^{(2)} \approx 0.65$, but various orientation parameters $r_N$ and $r_\alpha$. The time-evolution is computed following the two approaches of the previous “first” method (lines with points for the micro-macro approach, and continuous lines for Jeffery’s approach). A third approach of the time-evolution of the tested microstructures is also used (continuous lines with triangles): it consists in adopting the Jeffery’s approach to compute the orientation rate $\dot{\theta}_b^{[0]}$ of each fibre and incorporating the translation field $\epsilon v_b^{[1]}$ of each fibre in Eq. (32): this third approach therefore adopts a non-affine velocity field.

The graphs plotted in these figures conjure up the following comments:

* Whatever the initial fibrous microstructures and the method adopted to compute the evolution of fibre orientation, fibrous networks align in the $e_1$-direction during the flow: $\tilde{A}_{11}^{(2)}$ increases and $\tilde{A}_{12}^{(2)}$ remains very close to 0.

* As already pointed out in the previous section, when assuming randomly distributed fibres with Gaussian-based orientations (“fit” results), predictions given by Jeffery’s equation or by the micro-macro approach are very close, except for fibrous networks that are highly oriented in the $e_1$-direction (Fig. 9(a3) and (b3)). On the contrary, predictions given with the diffusion coefficient (PT “fit” results) systematically deviate from the two previous types of simulation: the fibre orientation is slowed down for initially random or aligned microstructures (Fig. 9(a2,b2,a3,b3)), but it is enhanced when fibres are initially aligned in the $e_2$-direction (Fig. 9(a1,b1)).

* Whatever the initial fibre orientation, predictions given by the first calculation method (micro-macro or Jeffery) always differ from those given by the second method (“fit” results, micro-macro or Jeffery). In order to analyze more closely such a discrepancy, we plotted in Fig. 10 the average orientation distribution functions that were obtained at the end of the simulations with the first method (micro-macro or Jeffery). We also plotted in the same figure the orientation distribution functions of REV’s having the same second-order orientation tensors but with Gaussian-based orientation distributions, i.e., similar to those that were used with the second method. As evident from this figure, the fibre orientation within flow-induced microstructures (method 1) severely differs from the Gaussian type orientation distribution: they may exhibit more than one preferential orientation peak (see Fig. 10(a)), and the corresponding peaks are sharper. It has to be noticed that this type of fibre orientation distributions have a shape quite close to the crenellated orientation distributions tested in the previous section. It has been observed that developing sharp orientation peaks in the fibre orientation distribution has for effect to decrease the orientation rate $\tilde{A}_{11}^{(2)}$. This observation may explain the strong decrease of the orientation rate $\tilde{A}_{11}^{(2)}$ observed in Fig. 9(b1-b3).

* For the investigated flow conditions, Fig. 9 also shows that the Jeffery’s equation and the micro-macro approaches for flow-

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

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

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Average orientation distributions for 20 different cells. Comparison of time-evolution calculations and initial Gaussian-based distributions with the same average orientation: (a) initial $\tilde{A}_{11}^{(2)} \approx 0.25$, (b) initial $\tilde{A}_{11}^{(2)} \approx 0.50$, (c) initial $\tilde{A}_{11}^{(2)} \approx 0.75$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Comparison of time-evolution calculations using the micro-macro and the Jeffery approaches, as well as the hybrid and natural closure approximations for 1D elongational flows in $P$ along $e_1$: (a) $\bar{A}_{11}^{(2)}$ with respect to time, (b) $\bar{A}_{11}^{(2)}$ with respect to $\bar{A}_{11}^{(2)}$, for initial random orientation ($\bar{A}_{11}^{(2)}(t=0) \cong 0.50$).</figcaption></figure>

induced microstructures give, at first approximation, nearly the same results concerning the second-order orientation tensor during the simulations. This remark is strengthened by the graphs of Fig. 10: the corresponding orientation distribution functions are rather close. The observed differences, which are very weak, are much pronounced for initially oriented microstructures (Figs. 9(b1) and 10(a)). They may be ascribed to a loss of the initial random positions of fibres when using the micro-macro approach, since in this case $\epsilon v_b^{[1]} \neq 0$.

* Fig. 11 shows that the natural closure approximation gives results very close to the Jeffery predictions. This was expected as this closure approximation was particularly developed to be exact when starting with an isotropic orientation, and if every fibre rotated according to Jeffery's equation. It must be pointed out that this result holds only for this particular case of fibre orientation and is not a general rule. On the contrary, the hybrid closure results are not in accordance with the evolution predicted

by the micro-macro or Jeffery-type approaches, the orientation state and the evolution rates are overestimated.

* Lastly, the continuous lines plotted in the graph of Fig. 12 prove that Jeffery's approach proposes reasonable predictions, at first approximation, of fibre orientation within the considered concentrated bundle suspensions even for initially crenellated fibre orientation distributions. By looking at the continuous line with the triangles, also notice that the prediction can be improved if the third approach (cf. above) of time-evolution simulations is used, i.e., when adopting a non-affine velocity field.

## 6. Conclusion

A micro-macro approach, i.e., a discrete homogenization technique, is used to study the rheological behaviour of very concentrated suspensions seen as connected planar networks of fibre-bundles with Newtonian viscous interactions. The numerical tool, which was developed on the basis of this theoretical approach, enables to generate fibre-bundle networks having various orientation distributions: i.e., Gaussian-based, uniform and crenellated orientation distribution functions. These networks can be subjected to any kind of mechanical loadings. In this work, 1D elongational and pure shear flows were used. It is thus possible to analyze the influence of microstructural parameters such as the fibre-bundle content or the fibre-bundle orientation on the evolution of the orientation state.

Results obtained from *instantaneous tests* first show that there is no influence of the fibre-bundle volume fraction on the orientation evolution in the concentrated regime, which is studied here. On the contrary, for a given type of orientation distribution, the mean average orientation given by the second-order orientation tensor, has a strong impact on the orientation rate. Changing the type of orientation distributions has also a drastic effect. For 1D elongational flows, a rule seems to emerge from the set of results: the sharper the peak exhibited by the orientation distribution, the lower the rate of orientation evolution. In spite of the completely different physics of both approaches, it is also shown that Jeffery's predictions can be well compared to the micro-macro results even if some discrepancies are found for microstructures aligned along the flow direction. Closure approximations are rather deficient. They are only relevant for some particular cases of orientation states. In the same way, the approach using a hydrodynamic diffusion term is deficient as well, and has also the disadvantage of losing the linearity of the rate of

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Evolution of $\bar{A}_{11}^{(2)}$ with $\bar{A}_{11}^{(2)}$ for initially crenellated fibre distributions (1D elongational flows in $P$ along $e_1$).</figcaption></figure>
the second-order orientation tensor with respect to the strain rate. This is not predicted by the micro-macro or Jeffery's approaches for the current case of a local physics based on viscous interactions between fibres.

"Full" micro-macro *time-evolution calculations* were performed. They consist in updating the position and orientation of each fibre-bundle of the tested microstructures. Results predict an evolution towards microstructures having orientation distributions which may exhibit sharp orientation peaks. Using Jeffery's equation to perform similar calculations for the rotation of fibre-bundles gives results that correlate quite well with those of the micro-macro approach. This was shown to work, when adopting an affine displacement field for each fibre-bundle, and was improved by applying also a non-affine one deduced from the micro-macro calculations. The good comparison between the micro-macro and Jeffery theories was observed for all tested orientation distributions. Thus, the Jeffery's theory seems accurate provided that an accurate fibre-bundle orientation is used at each calculation step. The little differences between both approaches for very oriented microstructures are an interesting point that remains to be further studied.

Results of the various instantaneous tests can also be fitted to give equations for the rate of evolution of the orientation of bundles. Solving these equations gives access to another method to calculate the time-evolution of fibre-bundle orientation. When adopting this method, it appears that predictions of approaches using closure approximations or an hydrodynamic diffusion term always largely differ from the micro-macro one. On the contrary, Jeffery's predictions appear to correlate quite well with those of the micro-macro one. Note also the very special case of the natural closure approximation, which matches the Jeffery prediction when starting from an initially isotropic orientation state. This result was expected as the natural closure approximation was built to fulfil this situation, but this is not a general result. Furthermore, it seems illusory to find a closure approximation that covers all encountered fibre orientation conformations, so that other solving strategies should be used to compute the evolution of fibre orientation.

Results gained by the “full” and “fitted” *time-evolution calculations* always differ even when performing a comparison between the micro-macro predictions for both methods. This result is not so surprising regarding the development of complex orientation states, and consequently complex fibre orientation distributions that appear when performing full time-evolution calculations of the first method. These distributions largely differ from the considered case of Gaussian-based orientation distributions used by the “fitted” method of time-evolution tests, which therefore cannot be used as a predictive technique.

Finally, in view of the complexity of fibre orientation states that can develop during the flow in concentrated suspensions, and of the importance of the fibre orientation distribution on the orientation rate, the solving strategy based on the use of closure approximation appears to be rather limited. Recently, a new discretization technique [39], which enables to reduce the number of degrees of freedom involved in the evolution equation for the orientation distribution function of type (6), has emerged. Such a strategy should preferentially be followed. This latter equation needs to incorporate a correct description of the motion of individual fibres, e.g., adapted to the physics of the considered suspension. As shown here, Jeffery's theory gives quite accurate predictions for the special case of the Newtonian concentrated suspensions, which have been studied. Future efforts should also focus on cases of concentrated suspensions where the suspending fluid and the interaction between fibres are non-linear. This point has been only slightly addressed in [40,7] and remains to be deeply investigated.

## References

[1] S. Le Corre, L. Orgéas, D. Favier, A. Tourabi, A. Maazouz, C. Venet, Shear and compression behaviour of sheet molding compounds, Compos. Sci. Technol. 62 (4) (2002) 571-577.

[2] C. Servais, A. Luciani, J.-A.E. Måanson, Squeeze flow of concentrated long fibre suspensions: experiments and model, J. Non-Newtonian Fluid Mech. 104 (2002) 165-184.

[3] P. Dumont, L. Orgéas, S. Le Corre, D. Favier, Anisotropic viscous behaviour of sheet molding compounds (SMC) during compression molding, Int. J. Plasticity 19 (4) (2003) 625-646.

[4] S. Toll, J.-A.E. Måanson, Dynamics of a planar concentrated suspension with non-hydrodynamic interaction, J. Rheol. 38 (4) (1994) 985-997.

[5] C. Servais, A. Luciani, J.-A.E. Måanson, Fiber-fiber interaction in concentrated suspensions: dispersed fiber bundles, J. Rheol. 43 (4) (1999) 1005-1018.

[6] S. Le Corre, D. Caillerie, L. Orgéas, D. Favier, Behavior of a net of fibers linked by viscous interactions: theory and mechanical properties, J. Mech. Phys. Solids 52 (2004) 395-421.

[7] S. Le Corre, P. Dumont, L. Orgéas, D. Favier, Rheology of highly concentrated fiber suspensions, J. Rheol. 49 (5) (2005) 1029-1058.

[8] C.L. Tucker, S.G. Advani, Processing of short-fiber systems, in: Composite Materials Series, volume 10, Elsevier, Amsterdam, 1994.

[9] F. Dupret, V. Verleye, Modelling the flow of fiber suspensions in narrow gaps, in: Rheology Series, Elsevier, Amsterdam, 1999, pp. 1347-1398.

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

[11] G.L. Hand, A theory of dilute suspensions, Arch. Rational Mech. Anal. 7 (1961) 81-86.

[12] J.L. Ericksen, Anisotropic fluids, Arch. Rational Mech. Anal. 4 (1960) 231-237.

[13] F. Folgar, C.L. Tucker, C.C. Lee, Simulation of compression molding for fiber-reinforced thermosetting polymers, J. Eng. Ind. - Trans. ASME 106 (1984) 114-125.

[14] S.G. Advani, C.L. Tucker, The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. 3 (8) (1987) 751-784.

[15] S. Prager, Stress-strain relations in a suspension of dumbbells, Trans. Soc. Rheol. 1 (1957) 53-62.

[16] S.G. Advani, C.L. Tucker, A numerical simulation of short fiber orientation in compression molding, Polym. Compos. 11 (3) (1990) 164-173.

[17] I. Doghri, L. Tinel, Micromechanics of inelastic composites with misaligned inclusions: numerical treatment of orientation, Comput. Meth. Appl. Mech. Eng. 195 (2006) 1387-1406.

[18] G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech. 13 (1962) 33-46.

[19] E.J. Hinch, L.G. Leal, Constitutive equations in suspension mechanics. Part 2. Approximate forms for a suspension of rigid particles affected by brownian rotations, J. Fluid Mech. 76 (1) (1976) 187-208.

[20] M. Rahnama, D.L. Koch, E.S.G. Shaqfeh, The effect of hydrodynamic interactions on the orientation distribution in a fiber suspension subject to simple shear flow, Phys. Fluids 7 (3) (1995) 487-506.

[21] E.S.G. Shaqfeh, D.L. Koch, Orientational dispersion of fibers in extensional flows, Phys. Fluids A 2 (7) (1990) 1077-1093.

[22] D.L. Koch, E.S.G. Shaqfeh, The average rotation rate of a fiber in the linear flow of a semidilute suspension, Phys. Fluids A 2 (12) (1990) 2093-2102.

[23] D.L. Koch, A model for orientational diffusion in fiber suspensions, Phys. Fluids 7 (8) (1995) 2086-2088.

[24] M.P. Petrich, D.L. Koch, C. Cohen, An experimental determination of the stress-microstructure relationship in semi-concentrated fiber suspensions, J. Non-Newtonian Fluid Mech. 95 (2000) 101-133.

[25] C.L. Tucker, R.B. Dessenberger, Flow and Rheology in Polymer Composites Manufacturing, Elsevier Sc, 1994.

[26] Y. Yamane, Y. Kaneda, M. Doi, The effect of interaction of rodlike particles in semi-dilute suspensions under shear flow, J. Phys. Soc. Jpn. 64 (9) (1995) 3265-3274.

[27] N. Phan-Thien, X.J. Fan, R.I. Tanner, R. Zheng, Folgar-tucker constant for a fibre suspension in a newtonian fluid, J. Non-Newtonian Fluid Mech. 103 (2002) 251-260.

[28] X. Fan, N. Phan-Thien, R. Zheng, A direct simulation of fibre suspensions, J. Non-Newtonian Fluid Mech. 74 (1998) 113-135.

[29] T.-H. Le, P.J.J. Dumont, L. Orgéas, 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.

[30] P.J.J. Dumont, J.P. Vassal, L. Orgéas, V. Michaud, D. Favier, J.-A.E. Måanson, Processing, characterisation and rheology of transparent concentrated fibre-bundle suspensions, Rheol. Acta 46 (2007) 639-651.

[31] H. Tollenaere, D. Caillerie, Continuous modelling of lattice structures by homogenization, Adv. Eng. Soft. 29 (7) (1998) 699-705.

[32] G. Moreau, D. Caillerie, Continuum modelling of lattice structures in large displacement framework, Dev. Comput. Tech. Struct. Eng. 20 (1995) 53-70.

[33] G. Moreau, D. Caillerie, Continuum modelling of lattice structures in large displacement - applications to buckling analysis, Comput. Struct. 68 (1998) 181-189.

[34] M. Doi, S.F. Edwards, Dynamics of rod-like macromolecules in concentrated solution, J. Chem. Soc. Faraday Trans. II 74 (1) (1978) 560-570.

[35] S. Ranganathan, S.G. Advani, Fiber-fiber interactions in homogeneous flows of non-dilute suspensions, J. Rheol. 35 (1991) 1499-1522.

[36] S. Toll, Note: On the tube model for fiber suspensions, J. Rheol. 37 (1) (1993) 123-125.
[37] T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: statistical and numerical approach, Int. J. Solids Struct. 40 (2003) 3647-3679.

[38] S. Rolland du Roscoat, M. Decain, X. Thibault, C. Geindreau, J.-F. Bloch, Estimation of microstructural properties from synchrotron X-ray microtomography and determination of the rev in paper materials, Acta Mat. 55 (2007) 2841-2850.

[39] A. Ammar, D. Ryckelynck, F. Chinesta, R. Keunings, On the reduction of kinetic theory models related to finitely extensible dumbbells, J. Non-Newtonian Fluid Mech. 134 (2006) 136-147.

[40] P. Dumont, S. Le Corre, L. Orgéas, and D. Favier Predicting the flow induced fibre orientation of highly concentrated suspensions using a micromacro modelling. In: V Bruscato (Ed.), 6th ESAFORM Conf. on Material Forming, pages 831-834, Salerno, Nuova Ipsa, 2003.