# On the Steady Regime Assumptions in Flows Involving Short Fiber Suspensions

J. Azaiez, K. Chiba, Francisco Chinesta, Amine Ammar

► To cite this version:

J. Azaiez, K. Chiba, Francisco Chinesta, Amine Ammar. On the Steady Regime Assumptions in Flows Involving Short Fiber Suspensions. International Journal of Forming Processes, 2006, 9 (4), pp.431-456. hal-01007157

HAL Id: hal-01007157

https://hal.science/hal-01007157

Submitted on 2 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# On the steady regime assumptions in flows involving short fiber suspensions

J. Azaiez<fnref n="1" />, K. Chiba<fnref n="2" />, F. Chinesta<fnref n="3" />, A. Ammar<fnref n="4" />

ABSTRACT. The present study addresses the issue of the steady regime assumption in flows involving short fibers usually introduced for the analysis of rheometric flows or in the simulation of some kind of forming processes as is the case of the extrusion. For this purpose we analyze the solution of the Fokker-Planck equation in simple shear flows for both two-dimensional and three-dimensional systems. Despite the simplicity of that flow, it brings the possibility of comparing exact and numerical solutions. Thus, exact analytical solutions and results of numerical simulations based on the finite element technique and the method of particles are presented. The effect of diffusion on the existence of steady solutions, and more particularly on the time required to reach that steady solution, is investigated.

KEYWORDS: Short fiber suspensions, Orientation distribution, Fokker-Planck equation.

<footnote n="1">¹ Department of Chemical and Petroleum Engineering
The University of Calgary,
2500 University Drive N.W., Calgary, Alberta, Canada T2N 1N4
azaiez@ucalgary.ca</footnote>

<footnote n="2">² Faculty of Education, Shiga University
2-5-1 Hiratsu, Otsu, Shiga, 520-0862 Japan
khiba@sue.shiga-u.ac.jp</footnote>

<footnote n="3">³ Laboratoire de Mécanique des Systèmes et des Procédés
UMRCNRS-ENSAM-ESEM
151 Boulevard de l'Hôpital, F-75013 Paris, France
francisco.chinesta@paris.ensam.fr</footnote>

<footnote n="4">⁴ Laboratoire de Rhéologie, INPG, UJF, CNRS (UMR 5520)
1301 rue de la piscine, BP 53 Domaine universitaire
38041 Grenoble Cedex 9, France
Amine.Ammar@ujf-grenoble.fr</footnote>
# 1. Introduction

Rigid particle suspensions are encountered in a wide variety of engineering applications related to materials processing. The orientation of the particles, which depends on the flow of the material during the manufacturing process, has a great influence on the physical properties of the final product (Pipes et al. 1982; Agari et al. 1991). A wealth of experimental and numerical studies has focused on modeling the flow of fiber suspensions in a variety of geometries (Rao et al. 1991; Keiller and Hinch 1991; Tang and Altan 1995; Ausias et al. 1996; Azaiez et al. 1997; Chiba and Nakamura 1998). These studies, whether adopting a suspending matrix which is Newtonian or non-Newtonian, showed conclusively the importance of determining exactly the orientation distribution of the particles in the flow. In particular, many numerical studies have attempted to determine the steady-state orientation distribution of fibers in simple as well as complex flow geometries due to its importance in many applications. Indeed, most rheometric devices that are commonly used to characterize the rheological behavior of fiber suspensions are based on the implicit assumption that the fiber suspension reaches a steady state during the operation time. Furthermore, quite few models of processes used in the manufacturing of composite materials such as extrusion focus on the steady-state solution that governs the properties of the processed material. In these systems and others, it is tacitly assumed that the steady-state orientation exists, that it is reached and perfectly established during the process, which serves as the main basis for the conclusions of the analysis and measurements.

The objective of the present study is to address two major issues. The first one deals with the existence of a steady-state orientation distribution in a simple shear flow for both two-dimensional and three-dimensional orientation distributions. In particular the effect of particle diffusion on the existence of such solution, as a long time solution, will be investigated. The second issue is related to the characteristic time required to reach this steady solution (when it can be reached as the long time solution of the associated evolution problem) with respect to the process time.

The present paper is organized as follows. The next section presents the mathematical model used to describe the flow of the suspension and to determine the fibers orientation distribution. In the third section, exact analytical solutions of the steady state orientation distribution in the presence as well as in the absence of orientation diffusion are obtained. Results of numerical simulations based on the finite element method as well as on particle techniques are presented in section four and compared with the analytical solutions of section three. An analysis and discussion of the results are presented in the last section.
## 2. Mechanical Modeling

The governing equations to describe the Stokes flow of a fiber suspension are the continuity and Cauchy momentum equations:

$$ \nabla \cdot \vec{v} = 0 \qquad [1] $$

$$ \nabla \cdot (-\mathbf{p}\mathbf{I} + \boldsymbol{\tau}) = 0 \qquad [2] $$

In the above equations, $\vec{v}$ is the velocity vector and $\boldsymbol{\tau}=\boldsymbol{\tau}^n+\boldsymbol{\tau}^f$ is the extra-stress tensor accounting for the contribution of the Newtonian matrix $\boldsymbol{\tau}^n$ and that of the fiber $\boldsymbol{\tau}^f$. The fiber contribution to the total stress is in general a function of the fiber properties (shape and dimensions), the fiber orientation distribution and the suspending fluid properties. It will be assumed here that the fibers are rigid particles with uniform length $\mathbf{I}$ and diameter $\mathbf{d}$, that they are homogeneously distributed in the suspension and have the same density as the suspending fluid.

The orientation of a single fiber in space can be represented by a unit vector $\vec{p}$ along the fiber axis or equivalently two angles $\theta$ and $\phi$ as shown in Figure 1. All existing models used to describe the orientation of fibers are based on Jeffery's equation for the motion of a single rigid ellipsoidal particle suspended in a Newtonian Stokes flow (Jeffery 1922):

$$ \dot{\vec{p}} = \Omega \vec{p} + k \mathbf{D} \vec{p} - k(\vec{p}^T \mathbf{D} \vec{p}) \vec{p} \qquad [3] $$

where $\Omega$ is the vorticity tensor, $\mathbf{D}$ the rate-of-strain tensor and $k=(r_e^2-1)/(r_e^2+1)$ is a parameter related to the aspect ratio of the fiber defined as the ratio of its length $\mathbf{I}$ and diameter $\mathbf{d}$: $r_e=\mathbf{I}/\mathbf{d}$. In most practical applications, a large number of fibers are suspended in the flow, and the best approach to determine the orientation distribution is based on a statistical description using an orientation distribution function $\Gamma(\vec{p})$, where $\Gamma(\vec{p})d\vec{p}$ is the probability that a fiber be oriented between $\vec{p}$ and $\vec{p}+d\vec{p}$. Owing to the complexity and the prohibitively high computational cost associated with the use of the statistical description particularly in complex flows, most existing studies have adopted an average-based approach using moments of the orientation distribution function $\Gamma(\vec{p})$. In particular, a second- and a fourth-order orientation tensor in terms of moments of the probability distribution function $\Gamma(\vec{p})$ are defined as (Advani and Tucker 1987):
$$ \mathbf{a}_2 = \langle \bar{\mathbf{p}}\bar{\mathbf{p}} \rangle = \int \bar{\mathbf{p}} \otimes \bar{\mathbf{p}} \Gamma(\bar{\mathbf{p}}) d\bar{\mathbf{p}} \qquad [4] $$

$$ \mathbf{a}_4 = \langle \bar{\mathbf{p}}\bar{\mathbf{p}}\bar{\mathbf{p}}\bar{\mathbf{p}} \rangle = \int \bar{\mathbf{p}} \otimes \bar{\mathbf{p}} \otimes \bar{\mathbf{p}} \otimes \bar{\mathbf{p}} \Gamma(\bar{\mathbf{p}}) d\bar{\mathbf{p}} \qquad [5] $$

<figure><img src="image_1.png" /><figcaption>Figure 1. Orientation of a single rigid fiber in space</figcaption></figure>

Evolution equations can be derived for the orientation tensors, and the numerical solutions of such equations are less computationally expensive than solving for the orientation distribution. Even though it is practical, the average-based approach has two major drawbacks. First, the orientation tensors are terms of an infinite series expansion of the distribution function $\Gamma(\bar{\mathbf{p}})$, and therefore it is not possible with only a few orientation tensors that contain only average-based information to recover a full and accurate description of the exact orientation distribution function. Second, the evolution equation for the second-order orientation tensor involves the fourth-order orientation tensor. In fact the evolution equation for the n-th order orientation tensor involves the (n+2)th order orientation tensor, and therefore an approximation that allows to express the (n+2)th order tensor in terms of lower-order tensors is necessary to solve the problem in a closed form. Different types of closure approximations particularly for the fourth-order orientation tensor $\mathbf{a}_4$ in terms of the second-order orientation tensor $\mathbf{a}_2$ have been proposed in the literature. Among the most commonly used closure approximations one may mention the quadratic closure (Doi and Edwards 1998), the linear closure (Hand 1962), the HL closure (Hinch and Leal 1975), the hybrid closure (Advani and Tucker 1990), the natural closure (Verleye and Dupret 1993), the orthotropic closure (Cintra and Tucker 1995) and the IBOF closure (Chuang and Kwon 1997). However existing studies show that the predictions of the average-based models can be strongly dependent on the type of closure approximations used, and so far none of the existing approximation is capable of providing an exact formulation valid for any type of flow.
In order to avoid the limitations of the average-based formulation, we will adopt in the present study the exact description based on the orientation distribution function $\Gamma(\bar{\mathbf{p}})$ that satisfies the Fokker-Planck equation:

$$ \frac{\partial \Gamma}{\partial t} + \nabla_{\bar{\mathbf{p}}} \cdot [\dot{\bar{\mathbf{p}}} \Gamma - \mathbf{D}_r \nabla_{\bar{\mathbf{p}}} \Gamma] = 0 \quad [6] $$

This is essentially a conservation equation in the orientation space and has a drift term $\nabla_{\bar{\mathbf{p}}} \cdot (\dot{\bar{\mathbf{p}}} \Gamma)$ and a diffusion term $\nabla_{\bar{\mathbf{p}}} \cdot (\mathbf{D}_r \nabla_{\bar{\mathbf{p}}} \Gamma)$. The orientation distribution function must satisfy the periodicity condition $\Gamma(\bar{\mathbf{p}}) = \Gamma(-\bar{\mathbf{p}})$ as well as the normality condition $\int \Gamma(\bar{\mathbf{p}}) d\bar{\mathbf{p}} = 1$. In equation [6], $\mathbf{D}_r$ is a diffusivity tensor accounting for fiber interactions. In a typical flow of fiber suspensions, one may have to deal with many types of inter-particle interactions that can influence the physics of the flow. These interactions can be associated with hydrodynamic or colloidal forces, excluded volume and friction or direct particle-particle interactions. In the present study we adopt the formulation proposed by Folgar and Tucker (1984) in which the diffusivity tensor $\mathbf{D}_r$ is assumed to be isotropic and is expressed as $\mathbf{D}_r = \mathbf{D}_r \mathbf{I} = \mathbf{C}_1 \gamma \mathbf{I}$, where $\gamma$ is the magnitude of the equivalent rate-of-strain and $\mathbf{C}_1$ is a scalar interaction coefficient. The coefficient $\mathbf{C}_1$, known as the Folgar-Tucker constant, is a measure of the intensity of fiber interactions in the suspension which are assumed to have an overall randomizing effect resulting in an isotropic interaction diffusivity tensor. It should be mentioned that even though some studies have attempted to include the anisotropic nature of the fiber-fiber interactions (Koch 1995; Fan et al. 1998), however such non-isotropic interactions will not be examined here.

### 3. Analytical solutions for a shear flow

We will first examine the steady-state solution of the Fokker-Planck equation in a simple steady shear flow, $\bar{\mathbf{u}} = (G_y, 0, 0)$. Assuming that such solution exists, the corresponding orientation distribution function for an isotropic diffusivity is governed by:

$$ \nabla_{\bar{\mathbf{p}}} \cdot (\dot{\bar{\mathbf{p}}} \Gamma - \mathbf{D}_r \nabla_{\bar{\mathbf{p}}} \Gamma) = 0 \quad [7] $$

In what follows, we discuss the solution of equation [7] for both two-dimensional and three-dimensional orientation distributions.

#### 3.1. Two-dimensional orientation distribution

For a two-dimensional orientation distribution, equation [7] reads:
$$ \frac{\partial(\omega\Gamma)}{\partial\varphi} - D_r \frac{\partial^2\Gamma}{\partial\varphi^2} = 0 \qquad [8] $$

where $\bar{p} = (\cos\varphi, \sin\varphi)$ and $\omega = d\varphi/dt$. In what follows, the two cases where the diffusivity is either zero or not will be examined.

### Case 1: $Dr=0$

In the case of zero diffusivity, the steady-state solution must satisfy the following equation:

$$ \Gamma(\varphi) = \frac{A_1}{\omega} \qquad [9] $$

where $A_1$ is a constant and $\omega = d\varphi/dt$ is obtained from equation [3] as:

$$ \omega = -G[1 + k - 2k\cos^2\varphi]/2 \qquad [10] $$

whose solution is given by:

$$ \tan(\varphi) = \sqrt{\frac{1-k}{1+k}} \tan[\sqrt{1-k^2}(\Upsilon - Gt/2)] \qquad [11] $$

the constant of integration $\Upsilon$ being: $\Upsilon = \frac{1}{\sqrt{1-k^2}} \tan^{-1} \left( \sqrt{\frac{1+k}{1-k}} \tan(\varphi_0) \right)$.

The exact expression of $A_1$ is obtained from the normalization condition $\int_0^{2\pi} \Gamma d\varphi = 1$, leading to:

$$ \Gamma(\varphi) = \frac{\sqrt{1-k^2}}{2\pi[1+k-2k\cos^2\varphi]} \qquad [12] $$

A question may arise as to whether this steady state orientation distribution can be actually obtained as the long-time solution of the time-dependent problem. Indeed, the most common approach in numerical modeling studies consists of treating the problem as transient and stepping it in time using standard algorithms. The steady state orientation distribution is assumed to be the long-time solution. However as we shall see, such an assumption can be incorrect and does not always ensure reaching the steady state distribution.
In the case of zero diffusion, the solution of the Fokker-Planck equation along the characteristics leads to:

$$ \frac{dt}{1} = \frac{d\varphi}{\omega} = -\frac{d\Gamma}{\Gamma d\omega/d\varphi} \qquad [13] $$

which implies that:

$$ \frac{d\omega}{\omega} = -\frac{d\Gamma}{\Gamma} \qquad [14] $$

or:

$$ \Gamma(\varphi) = \Gamma_o(\varphi_o) \frac{\omega_o}{\varphi} = \Gamma_o(\varphi_o) \frac{1+k-2k\cos^2\varphi_o}{1+k-2k\cos^2\varphi} \qquad [15] $$

The above result can be also obtained using physical arguments stating that, for zero diffusion, the probability of finding a fiber oriented between $\varphi$ and $\varphi+d\varphi$ is the same as that of finding the same fiber at an earlier orientation say at t=0, between $\varphi_o$ and $\varphi_o+d\varphi_o$, i.e.

$$ \Gamma(\varphi)d\varphi = \Gamma_o(\varphi_o)d\varphi_o \Leftrightarrow \Gamma(\varphi) = \Gamma_o(\varphi_o) \varphi_o \qquad [16] $$

Using equations [11] and [15] and the fact that $\Gamma_o(\varphi_o)$ is the orientation distribution function at some initial time t=0, one gets:

$$ \Gamma(\varphi, t) = \frac{(1-k^2)\Gamma_o(\varphi_o)}{1-k^2 \cos 2\tau + k(\cos 2\tau - 1)\cos 2\varphi + k\sqrt{1-k^2} \sin 2\tau \sin 2\varphi} \qquad [17] $$

where $\tau = Gt\sqrt{1-k^2}/2$.

It is easy to show that the above expression of $\Gamma(\varphi)$ satisfies the Fokker-Planck equation with zero diffusion. It can be seen from Eq. [17] or Eq. [15] that the steady state distribution given by Eq. [12] can be obtained as the long time solution of the transient problem, if and only if the initial orientation distribution (at t=0) is the steady state solution given by Eq. [12]. In particular an initial random orientation distribution ($\Gamma_o(\varphi_o)=1/2\pi$) will result in a time-oscillating orientation distribution with a period $T = 2\pi/G\sqrt{1-k^2}$ that never reaches a steady state.

### Case 2: $Dr \neq 0$

In the presence of diffusion, the solution of equation [8] is:
$$ \Gamma(\varphi) = e^{\int_{0}^{\varphi} \frac{\omega}{D_r} d\varphi'} \left[ A_0 + \frac{A_1}{D_r} \int_{0}^{\varphi} e^{-\int_{0}^{\varphi} \frac{\omega}{D_r} d\varphi'} d\varphi' \right] \quad [18] $$

where $A_0$ and $A_1$ are constants. Using equation [10], and requiring the solution to be normalized and periodic with period $\pi$, one gets the following expression of $\Gamma$:

$$ \Gamma(\varphi) = \frac{e^{Z(\varphi)}}{F_1} \left[ \int_{0}^{\varphi} e^{-Z(x)} dx + F_0 \right] \quad [19] $$

where:

$$ Z(x) = \frac{G(2x - k \sin 2x)}{4D_r} \quad [20] $$

$$ F_0 = \frac{e^{z(x)}}{1 - e^{z(x)}} \int_0^x e^{-z(x)} dx \quad [21] $$

$$ F_1 = \int_0^{2\pi} e^{Z(\varphi)} \left( \int_0^{\varphi} e^{-Z(x)} dx + F_0 \right) d\varphi \quad [22] $$

Note that the above solution tends towards the random distribution ($\Gamma=1/2\pi$) for an infinite diffusivity ($D_r \to \infty$).

It is worth investigating whether this steady solution can be obtained from solving the transient problem. It should be first mentioned that, from the parabolic character of the Fokker-Planck equation, it can be expected that if the diffusivity is non zero, the transient problem will converge towards a steady state solution, being the main question the characteristic time for reaching that steady solution. This question will be addressed numerically in section 4. Thus, for small but non-zero diffusivity, it is not clear if this steady state orientation distribution can be reached from the solution of the transient problem, at least in a reasonable time. In order to address this issue, the orientation distribution $\Gamma$ is expressed as a series expansion in terms of the diffusion coefficient:

$$ \Gamma(\varphi, t) = \Gamma_0(\varphi, t) + D_r \Gamma_1(\varphi, t) + D_r^2 \Gamma_2(\varphi, t) + \quad [23] $$

Introducing the above series expansion into the general Fokker-Planck equation, and identifying the different orders of $D_r$, it results:
$$ \left\{ \begin{array}{l} \text{Order 0} \Rightarrow \frac{\partial \Gamma_0}{\partial t} + \frac{\partial(\omega \Gamma_0)}{\partial \varphi} = 0 \\[1em] \text{Order 1} \Rightarrow \frac{\partial \Gamma_1}{\partial t} + \frac{\partial(\omega \Gamma_1)}{\partial \varphi} = \frac{\partial^2 \Gamma_0}{\partial \varphi^2} \\[1em] \vdots \\[1em] \end{array} \right. \qquad [24] $$

where $\Gamma_0(\varphi, t)$ is the solution of the Fokker-Planck equation when diffusion vanishes and it is given by Eq. [12]. The first order term $\Gamma_1(\varphi, t)$ satisfies:

$$ \frac{\partial \Gamma_1}{\partial t} + \omega \frac{\partial \Gamma_1}{\partial \varphi} = - \frac{\partial \omega}{\partial \varphi} \Gamma_1 + \frac{\partial^2 \Gamma_0}{\partial \varphi^2} \qquad [25] $$

Using the method of characteristics leads to:

$$ \frac{dt}{1} = \frac{d\varphi}{\omega} = \frac{d\Gamma_1}{-\Gamma_1 \frac{d\omega}{d\varphi} + \frac{\partial^2 \Gamma_0}{\partial \varphi^2}} \Rightarrow \left\{ \begin{array}{l} \frac{d\varphi}{dt} = \omega \\ \frac{d\varphi}{\omega} = \frac{d\Gamma_1}{-\Gamma_1 \frac{d\omega}{d\varphi} + \frac{d^2 \Gamma_0}{d\varphi^2}} \end{array} \right. \qquad [26] $$

The problem is then reduced to a first order linear ordinary differential equation:

$$ \frac{d\Gamma_1}{d\varphi} + \Gamma_1 \frac{\frac{d\omega}{d\varphi}}{\omega} = \frac{\frac{d^2 \Gamma_0}{d\varphi^2}}{\omega} \qquad [27] $$

whose solution is:

$$ \Gamma_1(\varphi) = \frac{c}{\omega} + \frac{\frac{d\Gamma_0}{d\varphi}}{\omega} \qquad [28] $$

where c is a constant of integration that can depend on $\varphi_0$ and $\varphi = \varphi(\varphi_0, t)$ is the solution of the first relation in Eq. [26] and is given by Eq. [11]. In fact $\varphi$ is the angle defining the orientation of a hypothetical fiber that was aligned in the direction $\varphi_0$ at the initial time $t=0$. Note that $\varphi_0$ is the orientation associated with $\varphi$ at the initial time. The constant of integration is determined for each trajectory from Eq. [28], as:
$$c(\varphi_0) = \omega(\varphi_0)\Gamma_1(\varphi_0) - \left.\frac{d\Gamma_0}{d\varphi}\right|_{\varphi_0} \qquad [29]$$

Thus, the expression of the first-order term in the series expansion is:

$$\Gamma_1(\varphi) = \frac{\omega(\varphi_0)\Gamma_1(\varphi_0) - \left.\frac{d\Gamma_0}{d\varphi}\right|_{\varphi_0}}{\omega(\varphi)} + \frac{\left.\frac{d\Gamma_0}{d\varphi}\right|_{\varphi_0}}{\omega(\varphi)} \qquad [30]$$

It can be noticed that if $\Gamma_1(\varphi)$ depends on $\varphi_0$ then the problem cannot reach a steady state. In fact, each angular coordinate $\varphi$ is reached for different trajectories related to different initial orientations $\varphi_0$ at different times. Thus, if the solution given by Eq. [30] depends on the considered trajectory ($\varphi_0$) then at position $\varphi$ the solution necessary evolves in time, and never reaches a steady state. The only possibility for reaching a steady state is that the integration constant $c$ is the same for all trajectories, i.e.:

$$\omega(\varphi_0)\Gamma_1(\varphi_0) - \left.\frac{d\Gamma_0}{d\varphi}\right|_{\varphi_0} = K \qquad [31]$$

where $K$ is a constant that can be determined using the normality condition:

$$\int_{0}^{2\pi} \Gamma_{1}(\varphi) d\varphi = \int_{0}^{2\pi} \left| \frac{K + \left.\frac{d\Gamma_{0}}{d\varphi}\right|}{\omega(\varphi)} \right| d\varphi = 0 \qquad [32]$$

where the fact that $\int_{0}^{2\pi} \Gamma_{0}(\varphi) d\varphi = 1$ has been taken into account. Thus, a steady state can be reached as the long time solution of the transient problem if and only if the following initial condition is considered (Eq. [31]):

$$\Gamma_1(\varphi_0) = \frac{K + \left.\frac{d\Gamma_0}{d\varphi}\right|_{\varphi_0}}{\omega(\varphi_0)} \qquad [33]$$

which can be shown to correspond to the steady state (as expected using the ergodic assumption).
In summary, the previous analysis based on a series expansion in terms of the diffusion coefficient showed that for weak diffusion, no steady state can be reached at long times using arbitrary initial distributions and therefore, for practical purposes the steady state orientation distribution cannot be obtained using the transient approach. It should be stressed that this result does not prove that the steady state cannot be reached from a transient resolution; it only indicates that it cannot be reached in a reasonable time, related to the validity of the series expansion considered.

## 3.2. Three-dimensional orientation distribution

### Case 1: Dr=0

For zero diffusivity, there are in general infinite steady solutions. Indeed, a steady solution can be defined on each Jeffery orbit (using the previous reasoning used in the 2D case) as soon as the probability distribution is defined with respect to the Jeffery orbits. This result is in perfect agreement with the conclusions of Leal and Hinch (1971). However, these steady solutions cannot be reached by solving the evolution problem from an arbitrary initial condition as reported by Manhart (2003).

### Case 2: Dr≠0

In the presence of very small but non zero diffusion, Leal and Hinch (1971) showed that there is a steady state solution, which when expressed in terms of the angular coordinates reads:

$$ \Gamma(\theta, \varphi) = \frac{J(\theta, \varphi) r^2}{(r^2 \sin^2 \varphi + \sin^2 \theta \cos^2 \varphi)^{3/2}} \quad [34] $$

where:

$$ J(\theta, \varphi) = [H(r) + \frac{K(r)}{C^2} + \frac{M(r)}{C^4}] F(C^2, r) \quad [35] $$

$$ C^2 = \frac{r^2 \sin^2 \varphi + \sin^2 \theta \cos^2 \varphi}{r^2 \cos^2 \theta \cos^2 \varphi} \quad [36] $$

The detailed expression of the parameters H, K, M and the function F(C²) can be found in Leal and Hinch (1971). The question related to time for reaching such a one solution will be addressed numerically in section 4.3.
## 4. Numerical simulation for steady shear flows

In this section, we will attempt to confirm the conclusions of the previous theoretical analysis and assess the validity of using the transient approach to determine the steady state orientation distribution. To this end, numerical solutions of the Fokker-Planck equation are obtained for two-dimensional and three-dimensional orientation distributions in the transient regime, even if the choice of the simulation parameters is far from the ones usually encountered in polymer processes.

### 4.1. Two-dimensional orientation distribution with diffusion effects

The advection-diffusion Fokker-Planck equation in a transient regime is:

$$ \frac{\partial \Gamma}{\partial t} + \frac{\partial(\omega\Gamma)}{\partial \varphi} - D_r \frac{\partial^2 \Gamma}{\partial \varphi^2} = 0 \quad [37] $$

This equation is integrated using the finite element method with a Streamline Upwind Petrov-Galerkin scheme (SUPG) to account for advection terms in the context of both explicit and implicit time discretizations. The initial orientation distribution is assumed to be isotropic, i.e. $\Gamma(t=0) = 1/2\pi$. The shear rate has been set to

$G = 1$, and the fiber factor $k = 0.6$. The computing time step is $2 \times 10^{-5}$ and successive plots of the distribution function are presented for time increments of 0.2 for two values of the diffusion factor: $D_r = 0.05$ and $D_r = 0.5$.

Figure 2 depicting the variation of $\Gamma$ for weak diffusion ($D_r = 0.05$) clearly shows a time dependent orientation distribution. After a transition time, the solution oscillates into a bounded domain that evolves slowly with time. It is difficult to reach definite conclusions regarding the long time evolution of the solution due to the inevitable numerical diffusion introduced by the discretization technique. On the other hand, the results of the second simulation with a large diffusion ($D_r = 0.5$) clearly shows that the steady state can be reached (Figure 3) quickly. As we shall show later, the smaller the diffusion coefficient the larger the characteristic time required to reach the steady state.
<figure><img src="image_2.png" /><figcaption>Figure 2. Orientation distribution evolution for Dr = 0.05, k = 0.6 and G = 1</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3. Orientation distribution evolution for Dr = 0.5, k = 0.6 and G = 1</figcaption></figure>

To confirm the results just obtained and to ascertain that they are independent of
the numerical technique adopted, the same problem is solved using a radically
different approach where the advection-diffusion Fokker-Planck equation is solved by
the method of particles (Ammar and Chinesta 2002; Ammar et al. 2006), after writing
the Fokker-Planck equation as a purely advection equation using the following
formulation:
$$ \frac{\partial \Gamma}{\partial t} + \frac{\partial(\hat{\omega} \cdot \Gamma)}{\partial \varphi} = 0 \qquad [38] $$

with

$$ \hat{\omega} = \omega - D_r \frac{1}{\Gamma} \frac{\partial \Gamma}{\partial \varphi} \qquad [39] $$

Now, one can proceed by considering the equivalent fiber rotation velocity $\hat{\omega}$ given by Eq. [39] instead of the real rotation velocity that results from Jeffery's equation. In usual particle approximations one assumes $\Gamma^n(\varphi) = \sum_{i=1}^{i=N} \frac{1}{N} \delta(\varphi - \varphi_i^n)$ where $N$ is the number of particles considered (each one with “a mass” of $1/N$) and $\varphi_i^n$ their orientations at time $t_n$. In order to apply Eq. [39] one needs to evaluate the distribution function gradient, and then the derivative of the Dirac's distribution $d\delta(\varphi - \varphi_i^n)/d\varphi$ that cannot be defined in a proper way. Thus, a smoothed approximation is required, from which it results:

$$ \Gamma^n(\varphi) = \sum_{i=1}^{i=N} \frac{1}{N} \zeta_{\varepsilon}(\varphi - \varphi_i^n) \qquad [40] $$

where we consider the following Dirac delta smoothing:

$$ \zeta_{\varepsilon}(x) = \frac{\zeta\left(\frac{x}{\varepsilon}\right)}{\varepsilon}, \quad \zeta(x) = \frac{e^{-x^2}}{\sqrt{\pi}} \qquad [41] $$

which verifies the normality condition, and where $\varepsilon$ is the smoothing parameter. If an initial isotropic orientation distribution is considered, the evolution of this orientation distribution can be computed using the following algorithm:

* Consider $N$ fibers whose orientations are defined at time $t_0 = 0$ by:

$$ \varphi_i^0 = (i-1) \times \frac{2\pi}{N}, \quad i \in [1, \dots, N] $$

* For each time step:

* Compute the rotation velocity $\dot{\varphi}_i^n \equiv \hat{\omega}(\varphi_i^n)$ for each fiber using Eq. [39].

* Update the fibers orientation: $\varphi_i^{n+1} = \varphi_i^n + \varphi_i^n \Delta t$

* Compute the fiber orientation distribution: $\Gamma^{n+1}(\varphi) = \sum_{i=1}^{i=N} \frac{1}{N} \zeta_{\varepsilon}(\varphi - \varphi_i^{n+1})$
When this approach is implemented to solve the evolution problem for a fiber suspension characterized by $k = 0.6$, with $D_r = 0.5$, in a shear flow with $G = 1$, it was found that the computed solution reaches quickly a stable steady state. For a diffusion coefficient $D_r = 0.05$ no steady state was reached in a reasonable time. The solution oscillates again into the same bounded envelop that evolves slowly with time. These results are in perfect agreement with the ones obtained using the stabilized finite element method as well as with the conclusions of section 3.

We should finally note that there are other approaches to account for the diffusion effects. These may consist in transforming the Fokker-Planck equation into a stochastic equation as is usually done in molecular modeling (Keunings 2004) or in introducing the diffusion effects using a random walk technique as considered by Chinesta et al. (2003). These approaches have not been adopted in this study since for strong diffusion effects, the number of realizations required to account properly for these effects is too high, and often results in large simulation CPU time.

We will now attempt to determine how the characteristic time required to reach a steady state varies with the diffusion coefficient. Let's consider again the Fokker-Planck equation

$$ \frac{\partial \Gamma}{\partial t} + \nabla \cdot [\omega \Gamma - D_r \nabla \Gamma] = 0 \qquad [42] $$

Using the method of separation of variables, the solution may be sought in the form $\Gamma(\varphi, t) = T(t)A(\varphi)$ leading to:

$$ \frac{dT}{dt} A = - \frac{d(\omega A)}{d\varphi} T + D_r \frac{d^2 A}{d\varphi^2} T \qquad [43] $$

If $T(t)A(\varphi) \neq 0$, one can write:

$$ \frac{dT}{dt} = - \frac{d(\omega A)}{d\varphi} + D_r \frac{d^2 A}{d\varphi^2} = \lambda \qquad [44] $$

which implies:

$$ \left\{ \begin{array}{l} T = e^{\lambda t} \\ -\frac{d(\omega A)}{d\varphi} + D_r \frac{d^2 A}{d\varphi^2} = \lambda A \end{array} \right. \qquad [45] $$

Integrating Eq. [45b] results in:
$$ \int_{0}^{2\pi} D_r \frac{d^2 A}{d\phi^2} d\phi - \int_{0}^{2\pi} \frac{d(\omega A)}{d\phi} d\phi = \lambda \int_{0}^{2\pi} Ad\phi \quad [46] $$

Now, taking into account the periodicity of A and its derivatives, one gets:

$$ \lambda \int_{0}^{2\pi} Ad\phi = 0 \quad [47] $$

Due to the normality of the solution, one of the eigenfunctions has a non-zero integral and in consequence its associated eigenvalue must be zero. On the other hand the eigenfunctions related to non-zero eigenvalues have a zero integral. It is also possible to show that all other non-zero eigenvalues are complex. Indeed, multiplying Eq. [45-b] by $\omega A$ and $dA/d\phi$ and integrating between 0 and $2\pi$ results in the following two equations, respectively:

$$ \lambda \int \omega A^2 d\phi = D \int \omega A \frac{d^2 A}{d\phi^2} d\phi \quad [48] $$

$$ 0 = \int \frac{d(\omega A)}{d\phi} \frac{dA}{d\phi} d\phi \Rightarrow -\int \omega A \frac{d^2 A}{d\phi^2} d\phi = 0 \quad [49] $$

The above two equations show that $\lambda \int \omega A^2 d\phi = 0$, which for non-zero eigenvalues implies that $\int \omega A^2 d\phi = 0$. Noting that $\omega$ does not change sign in the interval $[0, 2\pi]$, the previous equality is possible only if the eigenfunction A is complex, and therefore the corresponding eigenvalue must be also complex. In order to describe the temporal evolution of the solution one must analyze numerically the behavior of the largest real-part of all eigenvalues.

A finite difference discretization is used to solve numerically the periodic eigenvalue ODE (Eq. [45-b]). The mesh size $h_\phi$ was chosen such that numerical instabilities are avoided. The resulting eigenproblem is cast in a matrix form $MA=\lambda A$ and solved to determine N couples eigenvalue-eigenvector where N is the number of nodes used in the discretization. Thus, the solution of Eq. [42] can be written as:

$$ \Gamma = \sum_{i=1}^{N} \beta_n e^{\lambda_n t} A_n(\varphi) \quad [50] $$

where the coefficients $\beta_n$ must be chosen to satisfy the initial condition. It was possible to verify numerically that the real parts of all the complex eigenvalues are negative. Moreover, it was easy to verify numerically that the eigenfunction related to the zero eigenvalue corresponds, after normalization, to the steady state solution given
by Eq. [19]. It was also found that the largest real part approaches zero as the diffusion coefficient becomes small (see table 1). Since the characteristic time of the transient process scales like the largest real part of the eigenvalues, this implies that the time required for reaching the steady state increases as the diffusion coefficient decreases.

<table><caption>Table 1. Evolution of the real part of the largest non-zero eigenvalue with the diffusion coefficient for a fiber suspension characterized by $k=0.6$.</caption><thead><tr><th>D<sub>r</sub></th><th>Largest non-zero eigenvalue real part</th><th>Characteristic time</th></tr></thead><tbody><tr><td>5</td><td>-4.99988</td><td>0.2000</td></tr><tr><td>0.5</td><td>-0.51546</td><td>1.9400</td></tr><tr><td>0.05</td><td>-0.07738</td><td>12.923</td></tr><tr><td>0.005</td><td>-0.00886</td><td>112.86</td></tr><tr><td>0.0005</td><td>-0.00092</td><td>1086.9</td></tr><tr><td>0.00005</td><td>-0.00009</td><td>11111</td></tr></tbody></table>

## 4.2. Two-dimensional orientation distribution neglecting diffusion effects

The flow of a suspension where the diffusion effects are neglected (the diffusion coefficient being zero is examined in this section to determine the impact of the non existence of a steady state orientation distribution in coupled (kinematics-fiber orientation) and uncoupled simulations. In the coupled case, the extra-stress tensor in Eq. [2] can be expressed as

$$ \tau = 2\eta\mathbf{D} + \mu\phi\mathbf{D}:\mathbf{a}_4 \quad [51] $$

where $\eta$ is the viscosity of the suspending Newtonian fluid, $\phi$ the volume fraction of fibers and $\mu$ a material constant.

When diffusion effects are neglected, particle techniques are the best numerical approach since they introduce small amount of numerical diffusion. This technique has been successfully used to model fiber orientation particularly in the work of Chiba and Nakamura (1998). It consists of assuming several fibers at certain points such that their orientation represents the initial known orientation distribution. Thus, for example in the case of an initial isotropic orientation state, one may consider N fibers located at each point with the following orientations:

$$ \varphi_i^0 = (i-1) \times \frac{2\pi}{N}, \quad i \in [1, \dots, N] \quad [52] $$
The orientation of each fiber is updated with time along its trajectory using Jeffery's equation:

➤ Compute the velocity at each fiber location: $\vec{v}(\vec{x}_i^{-n-1})$

➤ Move the center of gravity of each fiber:

$$ \vec{x}_i^{-n} = \vec{x}_i^{-n-1} + \vec{v}(\vec{x}_i^{-n-1}) \Delta t $$

➤ Compute the gradient of velocity at that point from the velocity field obtained in the previous step: $\text{Grad } v|_{\vec{x}_i^{-n-1}}$

➤ Compute the fiber rotation velocity from Jeffery's

$$ \text{equation: } \varphi_i = \varphi_i^{n-1} (\text{Grad } v|_{\vec{x}_k^{-n-1}}, k, \varphi_i^{n-1}) $$

➤ Update the fiber orientation:

$$ \varphi_i^n = \varphi_i^{n-1} + \varphi_i \Delta t $$

➤ Compute the “discrete” fiber distribution:

$$ \Gamma^n(\varphi) = \sum_{i=1}^{i=N} \frac{1}{N} \delta(\varphi - \varphi_i^n) $$

where $\delta(\varphi - \varphi_i^n)$ denotes Dirac's delta distribution. In the previously described algorithm, the diffusion effects have been neglected. For a particular orientation distribution, the second order orientation tensor can be then computed:

$$ \left\{ \begin{aligned} a_{11}^n &= \sum_{i=1}^{i=N} \cos^2 \varphi_i^n \frac{1}{N} \\ a_{12}^n &= \sum_{i=1}^{i=N} \cos \varphi_i^n \sin \varphi_k^n \frac{1}{N} \\ a_{22}^n &= \sum_{i=1}^{i=N} \sin^2 \varphi_i^n \frac{1}{N} \end{aligned} \right. \qquad [53] $$

The flow of a suspension characterized by $D_r = 0, k = 0.6$ and $\phi\mu/\eta = 10$, between two parallel plates separated by a distance $H$ is analyzed. We must notice that the scenario analyzed in this section has not a real sense because Eq. [51] is only applicable in the case of large fiber aspect ratios, that is not the case when one is assuming $k = 0.6$. The computational domain is given by $(x^*, y^*) \in [0,10] \times [0,1]$, where $x^* = x/H$ and $y^* = y/H$, with the flow taking place in the x-direction and $y^* = 0$ representing the flow symmetry axis. Based on symmetry considerations, only the lower half part of the channel is considered. A parabolic velocity profile and an isotropic fiber orientation distribution on the inflow boundary $x^* = 0$ is assumed. The number of particles was approximately of 2500.
In the decoupled case where the flow kinematics is not affected by the presence of fibers, all the streamlines are parallel to the x-axis. It is possible to integrate Jeffery's equation using the method of characteristics with a fourth order Runge-Kutta method along the flow streamlines. Figure 4 depicts the fiber orientation distribution at points $(x^*, y^*) = (0.5, 2/40)$ and $(x^*, y^*) = (0.7, 2/40)$. As expected, the fiber orientation distribution evolves periodically and is unable to reach a fully developed orientation regime. Since the period depends on the shear rate, the evolution of the fiber orientation distribution changes from one streamline to another, inducing large spatial heterogeneities.

<figure><img src="image_4.png" /><figcaption>Figure 4. Fiber orientation distribution in a flow between two parallel plates for</figcaption></figure>

$D_r = 0, k = 0.6$ and $y^* = 2/40$

<figure><img src="image_5.png" /><figcaption>Figure 5. Evolution of orientation ellipses along the streamlines defined by $y^*=1/40$, $2/40$ and $3/40$ at the channel inlet</figcaption></figure>

The spatial heterogeneities induced by the periodic evolution of the fiber orientation distribution along each streamline as well as the dependence of the period on the shear rate, and in consequence on the streamline, are expected to affect significantly the flow kinematics as one proceeds by coupling the kinematics of the
flow and the fiber orientation distribution. Figure 5 depicts the evolution of the orientation ellipses computed from the expression of the second order orientation tensor (Eq. [53]) in the coupled case. The direction and magnitude of the axes correspond with the principal directions and the eigenvalues of that tensor. The evolution of the orientation ellipses along the streamlines passing through y*=1/40, 2/40 and 3/40 at the channel inlet are presented. Note that the existence of a steady state in simple shear flows implies the existence of fully developed solutions. One possibility to circumvent this difficulty and reach a fully developed flow lies in the use of direct simulation considering a distribution of particles aspect ratio (Ausias, 2006).

### 4.3. Three-dimensional orientation distribution with diffusion effects

The Fokker-Plank equation is solved in the 3D case in a tessellation of the unit sphere surface, using the finite element method with the SUPG formulation as in the 2D case. We consider an initial isotropic fiber orientation distribution $\Gamma(t=0) = 1/4\pi$. The shear rate has been set to $G = 1$, and the fiber factor to $k = 0.6$. Two values of the diffusion coefficient are used $Dr_1 = 10^{-5}$ (for the first simulation) and $Dr_2 = 0.5$ (for the second one) to explore its effect on the long time solution in transient simulations.

The evolution of the fiber orientation distribution with $Dr_1 = 10^{-5}$ depicted in figure 6 on the unit surface proves that a steady state solution is not reached in a reasonable simulation time from the initial isotropic orientation distribution even if, still based on Leal and Hinch (1971) study, a steady state solution exists. In that figure the x-axis represents the flow direction being the xy-plane the shear plane. After a transition time, the solution oscillates with an amplitude that varies very slowly in time. To illustrate clearly this behavior, the evolution of the tensor $S$ related to the second moment of the distribution function is computed:

$$S = \int \bar{p} \otimes \bar{p} \Gamma(\bar{p}) d\bar{p} - I/3 \quad [54]$$

as well as the associated order parameter $s$:

$$s = \sqrt{\frac{3}{2} S : S} \quad [55]$$
<figure><img src="image_6.png" /><figcaption>Figure 6. 3D fiber orientation distribution (Dr = 10⁻⁵, k = 0.6, G = 1) at times t=1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12.</figcaption></figure>

Figures 7 and 8 show the evolution of $s$ and $S_{12}$ respectively, from which the oscillatory behavior of the orientation distribution is confirmed.
<figure><img src="image_7.png" /><figcaption>Figure 7. Order parameter evolution for $Dr = 10^{-5}$, $k = 0.6$ and $G = 1$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8. Evolution of $S_{12}$ for $Dr = 10^{-5}$, $k = 0.6$ and $G = 1$</figcaption></figure>

In the second simulation with $Dr = 0.5$, the steady state solution shown in the figure 9 is reached after a short transient period. This steady state solution does not depend on the initial distribution considered.
<figure><img src="image_9.png" /><figcaption>Figure 9. 3D steady state solution for Dr = 0.5, k = 0.6 and G = 1</figcaption></figure>

The above results indicate that for the three dimensional orientation, a steady state solution can be reached in a short time as long as the diffusion term is large enough. For weak diffusion, the steady state solution obtained by Leal and Hinch (1971) is not reached for a reasonable simulation time. It is however difficult to conclude about the behavior for extremely large times because of possible numerical diffusion introduced by numerical algorithms used to solve the problem, moreover when the simulation time becomes too large the inevitable numerical dissipation can induce non-physical behaviors.

5. Conclusion

Two major issues have been addressed in the present study. The first one deals with the existence of a steady-state orientation distribution in a simple shear flow for both two-dimensional and three-dimensional orientation distributions. In particular the effect of particle diffusion on the existence of the solution has been investigated. The second issue is related to the characteristic time required to reach one a such steady solution as the long time solution of a transient problem.

For the two-dimensional case, it was found that in the absence of diffusion, though a steady solution exists, it cannot be reached by solving the evolution problem from an arbitrary initial condition. The existence of a steady solution is also proved when the diffusion effects are included. In the case of very small diffusion effects, it was found that the characteristic time to reach the steady state increases as the diffusion
decreases, and therefore large computational times are required to get the steady state solution when the diffusion is very small.

In the case of a three-dimensional orientation, an infinite number of steady solutions exist when the diffusion is neglected but they cannot be obtained as the long time solution of a transient problem. When the diffusion is considered, numerical tests indicate that, as in the two-dimensional case, the steady solution can be reached in a reasonable time if an enough amount of diffusion is introduced.

## 6. References

Advani S.G., Tucker C.L., « The use of tensors to describe and predict fiber orientation in short fiber composites », *J. Rheol.*, vol. 31, 1987, p. 751-784.

Advani S.G., Tucker C.L. « Closure approximations for three-dimensional structure tensors », *J. Rheol.*, vol. 34, 1990, p. 367-386.

Agari Y., Ueda A., Nagai S., « Thermal conductivity of a polyethylene filled with disoriented short-cut carbon fibers », *J. Appl. Polym. Sci.*, vol. 43, 1991, p. 1117-1124.

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

Ammar A., Chinesta F., Ryckelynck D., « Deterministic Particle Approach of Multi-Bead-Spring Polymer Models », *European Journal of Computational Mechanics*, 2006, in press.

Ausias G., J. Jarrin J., Vincent M., « Optimization of the tube-extrusion die for short-fiber-filled polymers », *Composites Science and Technology*, vol. 56/7, 1996, p. 719-724.

Ausias G., Fan X.J., Tanner R., « Direct simulation for concentrated fiber suspension in transient and steady state shear flows », *Journal of Non-Newtonian Fluid Mechanics*, vol. 135, 2006, p. 46-57.

Azaiez J., Guénette R., Ait-Kadi A., « Investigation of the abrupt contraction flow of fiber suspensions in polymeric fluids », *J. Non-Newtonian Fluid Mech.*, vol. 73, 1997, p. 289-316.

Chiba K., Nakamura K., « Numerical solution of fiber suspension flow through a complex channel », *J. Non-Newtonian Fluid Mech.*, vol. 78, 1998, p. 167-185.

Chinesta F., Chaidron G., Poitou A., « On the solution of Fokker-Planck equations in steady recirculating flows involving short fiber suspensions », *Journal of Non-Newtonian Fluid Mech.*, vol. 113, 2003, p. 97-125.

Chuang D.H., Kwon T.H., « Improved model of orthotropic closure approximation for flow induced fiber orientation », *Polymer Comp.*, vol. 22, 2001, p. 636-649.
Cintra J.S., Tucker C.L. « Orthotropic closure approximation for flow-induced fibre orientation », **J. Rheol.**, vol. 39, 1995, p. 1095-1122.

Doi M., Edwards S.F. « Theory of Polymeric Dynamics », Oxford, 1998.

Fan X.J., Phan-Thien N., Zheng R., « A direct simulation of fibre suspensions », J. Non-Newtonian Fluid Mech., vol. 74, 1998, p. 113-135.

Folgar F.P., Tucker C.L., « Orientation behavior of fibers in concentrated suspensions », *J. Reinf. Plast. Comp.*, vol. 3, 1984, p. 98-119.

Hand G.L. « A theory of anisotropic fluids », **J. Fluid Mech.**, vol.13, 1962, p. 33-46.

Hinch E.J., Leal L.G., « Constitutive equations in suspension mechanics. Part 1. General formulation », **J. Fluid Mech.**, vol. 71, 1975, p. 481-495.

Jeffrey G.B., « The motion of ellipsoidal particles immersed in viscous fluid », *Proc. R. Soc. London*, vol. A 102, 1922, p. 161-179.

Keiller R.A., Hinch E.J., « Corner flow of a suspension of rigid rods », *J. Non-Newtonian Fluid Mech.*, vol. 40, 1991, p. 323-335.

Keunings R., in Rheology Reviews 2004, D.M. Binding and K. Walters (Eds.), British Society of Rheology, 2004, p. 67-98.

Koch D.L., « A model for orientational diffusion in fiber suspensions », **Phys. Fluids**, vol. 7, 1995, p. 2086-2088.

Leal L.G., Hinch E.J., « The effect of weak Brownian rotations on particles in shear flow » *J. Flui Mech.*, vol. 46, 1971, p. 685-703.

Manhart M., « Rheology of suspensions of rigid-rod like particles in turbulent channel flow », *Journal of Non-Newtonian Fluid Mech.*, vol. 12, 2003, p. 269-293.

Pipes R.B., Mc Cullough R.L., Taggart D.G., « Behavior of discontinuous fiber composites: Fiber orientation », *Polym. Compos.*, vol. 3, 1982, p. 34-39.

Poitou A., Chinesta F., Torres R., « Numerical simulation of the steady recirculating flows of fiber suspensions », *J. Non-Newtonian Fluid Mech.*, vol. 90, 2000, p. 65-80.

Rao B.N., Akbar S., Altan M.C., « Comparative study on the solution techniques for fiber orientation in two-dimensional converging and diverging flows », *J. Thermopl. Comp. Mat.*, vol. 4, 1991, p. 311-348.

Tang L., Altan M.C., « Entry flow of fiber suspensions in a straight channel », *J. Non-Newtonian Fluid Mech.*, vol. 56, 1995, p. 183-216.

Verleye V., Dupret F., Proc. ASME Winter Annual Meeting, New Orleans, LA, AMD-Vol. 175, p. 139-163, Developments in Non-Newtonian Flows, ASME November 28-December 3, 1993.