# State-of-the-Art on numerical simulation of fiber-reinforced Thermoplastic Forming Processes

Jalel Azaiez, Kunji Chiba, Francisco Chinesta, Arnaud Poitou

► To cite this version:

Jalel Azaiez, Kunji Chiba, Francisco Chinesta, Arnaud Poitou. State-of-the-Art on numerical simulation of fiber-reinforced Thermoplastic Forming Processes. Archives of Computational Methods in Engineering, 2002, 9 (2), pp.141-198. 10.1007/BF02736650. hal-01007319

HAL Id: hal-01007319

https://hal.science/hal-01007319

Submitted on 2 Feb 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# State-of-the-Art on Numerical Simulation of Fiber-Reinforced Thermoplastic Forming Processes

## J. Azaiez

Department of Chemical and Petroleum Engineering
The University of Calgary
2500 University Drive N.W.
Calgary, Alberta, Canada T2N 1N4
azaiez@ucalgary.ca

## K. Chiba

Faculty of Education, Shiga University
2-5-1 Hiratsu, Otsu, Shiga, 520-0862 Japan
kchiba@sue.shiga-u.ac.jp

## F. Chinesta

Laboratoire de Mécanique des Systèmes et des Procédés
ENSAM Paris, 151 Boulevard de l'Hôpital, F-75013 Paris, France
chinesta@paris.ensam.fr

## A. Poitou

LMT Cachan
61 Avenue du President Wilson
F-94235 Cachan cedex, France
poitou@lmt.ens-cachan.fr

Short fiber reinforced composites have gained increasing technological importance due to their versatility that lends them to a wide range of applications. These composites are useful because they include a reinforcing phase in which high tensile strengths can be reached, and a matrix that allows to hold the reinforcement and to transfer applied stress to it. It is a well-known fact that such materials can have excellent mechanical, thermal and electrical properties that make them widely used in industry. During the manufacture process, fibers adopt a preferential orientation that can vary significantly across the geometry. Once the suspension is cooled or cured to make a solid composite, the fiber orientation becomes a key feature of the final product since it affects the elastic modulus, the thermal and electrical conductivities, and the strength of the composite material. In this work we analyzed the state-of-the-art and the recent developments in the numerical modeling of short fiber suspensions involved in industrial flows.

# 1 INTRODUCTION

## 1.1 The Industrial Motivation

There are two main classes of composite materials depending on the field of their applications. In the first category one finds the composites dedicated to high diffusion parts in automobile industry, packing, house ware, etc. Their price exceeds rarely few dollars per kg. On the other hand, technical composites are widely used in sports equipment, and aerospace and naval industries due to their specific resistance in spite of their more expensive price. In most cases, technical composites consist of continuous fibers in a thermoset matrix, whereas high diffusion composites are very often made of short fibers reinforced thermoplastics (e.g. polypropylene, polyamide, etc.). In a volume average, the high diffusion composites represent the most essential part of the market. Many structural analysis studies concern high performance composite. On the contrary, it is understandable that many research studies
on forming processing are focused on high diffusion composites for which a process optimization is capital. This paper lies in this framework. Its aim is to establish the state of the art on the numerical simulation of short fibers reinforced thermoplastics (SFRT) flows, which are involved in short fiber composites forming processes.

There are many different forming processes for high diffusion composites. From a general point of view, every process used in thermoplastics transformation can be adapted in a more or less straightforward way to short fiber composites forming processes. However, among all these processes, injection molding and extrusion are probably the most usual. Moreover during these processes the material flows and induces a very large material deformation, which leads for composites, to a specific phenomenon: the fibers affect the flow and simultaneously, the flow induces a fiber orientation. These processes turn out to be a way to induce an anisotropy in the manufactured parts. This anisotropy is usually desired but has to be controlled in order to optimize the mechanical properties of the final product (fibers orientated locally in the stress directions).

**Injection molding process.** A typical injection molding machine is described in Figure 1. The process is sequential: the premixed material is introduced in a granular form into the barrel. The fusion arises because of the combined effects of electrical heating devices around the barrel and of internal dissipation effects due to the high viscosity of thermoplastics. The screw rotates and moves back to store the desired quantity of molten material at the head of the screw where the injection gate is located. Then, the rotation of the screw is stopped and the screw used as a piston, injects the material into the mold. The final part is then obtained after cooling under pressure. This technique is used to manufacture 3D parts.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Injection molding machine</figcaption></figure>

**Extrusion process.** Extrusion is used to obtain “one-dimensional” parts such as tubes, profiles and wires. It is basically the same process as injection molding with the exception that the process is continuous: the molten material is forced through a die and then generally cooled outside the extruder into a water bath (Figure 2). Changing the shape of the product cross section requires only to change the die geometry. Extruders work generally during many hours in a steady state regime.

**Induced anisotropy.** From a general point of view, processes concerning short fibers suspensions involve two main common features:

* A flowing device (mold, die, etc.), which orients the fibers.

* The presence of fibers affects the material behavior and thus the flow kinematics.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Extrusion machine</figcaption></figure>

*Flow orients the fibers.* To illustrate the first phenomenon from a physical point of view, we consider the injection of a disc. In our material, the fibers have a length around 0.1 mm and a diameter around 1µm. The fibers aspect ratio defined as the ratio of the length and diameter is equal to 100. The injection gate is located at the center of the disc on its surface. After the molding one cuts the part along a plane containing the disc axis and one of its diameter. If we examine the cross section with an optical microscope (see Figure 3), a clear heterogeneity of the morphology through the thickness can be detected. The fibers are found to be almost aligned with the flow near the upper and lower surfaces of the mold. On the contrary in the core region, near the middle plane of the disc, fibers are found to be mostly perpendicular to the flow. Between these two regions there is a transition zone where fibers seem to be oriented at random without any preferential alignment.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Fiber orientation distribution in an injected disc: experimental measurements (Meslin, 1997)</figcaption></figure>

These micrographs are an experimental evidence of the above-mentioned point of the effects of the flow kinematics on the fibers orientation. Before presenting an accurate modeling of this phenomenon, it would be interesting to give a simple semi-heuristic argument that supports these observations. Let's consider a simple shearing flow as shown in Figure 4 (fully-developed flow through a parallel plate channel). In this case, the fibers located out of the middle plane are subjected to a hydro-dynamical torque that induces a rotation of the fibers. If the fibers are assumed to have a quasi-infinite aspect ratio, they will rotate until they reach an equilibrium state in which they are fully aligned with the velocity field. The alignment characteristic time decreases when both (i) the flow rate decreases and (ii)
the fibers are located near the middle plane. When the fibers are very close to the flow symmetry zones, there is namely no velocity gradient and consequently, the hydro-dynamical torque vanishes. We can conclude that shear flows tend to align the fibers in the flow direction.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Orientation process in shear flow</figcaption></figure>

Lets consider now flows with an elongational component. A first example is a converging flow. An intuitive graphic argument shows that the torque due to the velocity gradient tends to align the fibers in the flow direction (Figure 5). On the contrary in a diverging flow, the torque aligns the fibers perpendicularly to the flow (Figure 6).

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Fiber orientation in convergent flow</figcaption></figure>

From these qualitative conclusions, the experimental results concerning the disc micrographs can be explained. In the zone near the walls the shear component is preponderant and the fibers tend to align in the flow direction. In the core region, the shear components almost vanish giving rise to a dominant extensional flow pattern similar to the one found in a divergent die (the norm of the velocity decreases along the trajectories): in this case the fibers align in the ortho-radial direction. In the transition region both the shear and elongational components act simultaneously giving rise to a non preferential fiber orientation.

*Effects of fibers orientation on the flow.* To show an experimental evidence of the influence of fibers on the flow pattern, lets consider a contraction flow at very low Reynolds number. Flow visualization with a planar sheet of laser light has have been achieved by Lipscomb et al. (1988). Figure 7 shows the case of a flow without fibers, where a typical small recirculation region can be seen. Figure 8 shows the case of the same material in which a very small amount of fibers has been added. It is very clear from this figure that the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Fiber orientation in divergent flows</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Newtonian fluid in a contraction flow (Lipscomb *et al.* (1988))</figcaption></figure>

fibers presence strongly affects the flow. A large recirculation region appears that cannot be explained in the framework of low Reynolds number Newtonian fluid mechanics. This phenomenon is somehow similar to the one encountered in viscoelastic fluids and can induce flow instabilities.

**Need for prediction and control of induced anisotropy.** Even if it is possible, as outline above to predict the fiber orientation in simple flows, in flows with industrial interest, this is not anymore possible: the complexity of the geometry as well as the coupling between the flow and fiber presence and orientation makes impossible an accurate prediction of the orientation. Without an accurate fiber orientation prediction it is not possible to evaluate the mechanical properties of the conformed part. One possibility to circumvent this difficulty consists in the experimental measurement of the fiber orientation on prototype pieces. However, such experimental trial-and-error approaches are known to be very expensive because they require not only the design but also the manufacture of several molds. Thus, the numerical simulation turns to be a very attractive solution from both theoretical and economical points of view.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Short fiber suspension in a contraction flow (Lipscomb *et al.* (1988))</figcaption></figure>

The aim of this paper lies in this framework of a need for a fully coupled modeling involving both the modeling of the materials constitutive behavior and the numerical modeling of the resulting equations.

With a forming process, one obtains usually an orientation field which is not optimal. The definition of an optimal criterion is in itself a rather complicated problem (limit analysis, fracture mechanics, resistance in fatigue or creep ...). However, assuming that we know the fiber orientation that we wish to obtain (stress directions), it can be possible to design the mold and the injection process, or the die geometry in extrusion processes, in order to fulfill the required orientation conditions. As an example, let's consider the case of an extruded tube (Figure 2). It can be interesting to orient the fibers on the circumferential direction in order to reinforce the tube to internal pressure loading conditions. If we design a die of Poiseuille's type, the flow is dominated by shear, which aligns the fibers on the tube axis direction (i.e. in the wrong direction from the mechanical point of view). However, using an axisymmetric divergent die the circumferential fiber orientation is improved.

## 1.2 State of the Art

### 1.2.1 *A general overview*

A mathematical modeling of a short fiber suspension flow involves mainly two steps. In the first step, one needs to derive a constitutive relation for the suspension. In the second one, it is necessary to develop accurate numerical schemes to treat the resulting equations.

Most models of anisotropic suspensions involve averaging procedures because it is not possible to follow the movement of each fiber in the suspension. These averaging procedures are of two kinds. The first one is a spatial averaging (often referred to as homogenization). It allows to derive macroscopic quantities as soon as the microscopic geometry is characterized. The second one is statistical: because the morphology is never fully characterized with exactitude, this method allows to derive statistical averaged macroscopic quantities from a statistical description of the morphology (Batchelor (1970-a), (1970-b), Hinch and Leal (1975) (1976), Meslin (1997)).

These homogenization procedures have been developed under different assumptions: the ambient fluid is Newtonian, there is not contact between fibers, fibers are rigid, inertia and
mass terms can be neglected, etc. Strictly speaking, these hypotheses are not realistic in real industrial situations where the matrix is visocelastic and the usually high volume fractions induces fiber contacts. However if one is not interested with a fine stress prediction it is admitted nowadays, that these oversimplified hypotheses allow for a reasonable prediction of the fiber orientation if the rheological parameters are identified experimentally.

The numerical modeling involves the coupling of an elliptic problem (anisotropic extension of Stokes equations) with a hyperbolic equation governing the fiber orientation evolution inside the suspension. Until now, there has been plane strain 2D, axisymmetric or 3D plate's type simulations in injection and extrusion processes. The 3D plate's type simulations are particularly interesting because many injection molding applications involve pieces which exhibit a dimension much smaller than the others. Several commercial codes already offer the possibility to solve the coupled problem for this type of geometries. Their numerical predictions are in good agreement with experimental results. However, some particular situations require real 3D calculations (small pieces or non axis-symmetric dies). Moreover simplified analysis leads to significant errors near the boundaries (fountain flow effect for example). Due to the hyperbolic nature of the orientation equation, slight errors are propagated in the domain by the flow. For the moment, there are no industrial simulation codes for these complex 3D flow geometries.

### 1.2.2 A more detailed overview: related works

Extensive numerical studies on the fiber suspension flows have been devoted to extrusion flow, flow around a sphere, squeezing flow, developing flow through a channel, radial flow between parallel disks and mold filling flow:

* Papanastasiou and Alexandrou (1987) investigated the isothermal extrusion of a non-dilute fiber suspension in a Newtonian solvent using the Dinh-Armstrong constitutive equation (Dinha and Armstrong (1984)). The fiber orientation was computed by solving the orientation distribution function along selected streamlines of the complex velocity field. Rosenberg et al. (1990) studied both planar extrusion and falling-ball rheometry of dilute fiber suspensions in a Newtonian fluid. They used the Transversely Isotropic Fluid (TIF) model first proposed by Ericksen (1960). The fiber orientation was obtained by numerically integrating the evolution equation of the second-order orientation tensor with a closure approximation along the streamlines.

* Phan-Thien and Graham simulated the squeezing flow between two circular disks of a model suspension (Phan-Thien and Graham (1990)) and the flow past a sphere (Phan-Thien and Graham (1991)). The suspension assumed to be a dilute and consisting of rigid, large aspect ratio spheroids, was adequately modeled by the TIF model.

* Altan et al. (1992) investigated the developing flow and fiber orientation of a suspension in a straight channel. They assumed that the orientation was planar and used the Dinh-Armstrong model. The orientation state was obtained by computing the evolution equation of the fourth-order orientation tensor with a sixth-order quadratic closure approximation. On the other hand, Chono and Makino (1995) used the Dinh-Armstrong model and the evolution equation of an equivalent strain tensor derived by Lipscomb (1986). Furthermore, Ausias et al. (1994) used the constitutive equation developed by Lipscomb et al. (1988) and the evolution equation of the fourth-order orientation tensor with a quadratic closure approximation. They computed an axisymmetric flow through an annular tube die by coupling with three-dimensional orientation of fibers.
* Ranganathan and Advani (1993) studied an axisymmetric diverging radial flow between disks of non-dilute fiber suspension. They coupled a modified form of the rheological equation given by Shaqfeh and Fredrickson (1990) with a model for the evolution of fiber orientation which takes into account the effect of fiber-fiber interactions on the fiber angular velocities. Altan and Rao (1995) developed the closed-form solution of the three-dimensional orientation field induced by a radially diverging, steady Newtonian flow between two parallel disks. Amhed and Alexandrou (1994) investigated the filling of a simple two-dimensional injection molding cavity. The Dinh-Armstrong model was used and the fiber orientation was computed by solving the orientation distribution function along the streamlines. Chung and Kwon (1995) predicted the transient behavior of fiber orientations together with a mold filling simulation of short-fiber suspension in arbitrary three-dimensional mold cavities. The Dinh-Armstrong model was incorporated into the Hele-Shaw equation and the evolution equation of the second-order orientation tensor was computed to obtain the orientation of fibers.

All of the above studies are concerned with non-recirculating flows and there has been relatively little work on fiber suspension flows including a recirculating flow as well as fiber orientation in a recirculating flow.

* Lipscomb et al. (1988) and Chiba et al. (1990) studied the circular entry flow of fiber suspensions in a Newtonian solvent. The TIF model was applied to suspensions with large aspect-ratio fibers. The results obtained using this model were compared with experimental data and the growth of the salient corner vortex was predicted with high accuracy over the semi-dilute fiber suspension regime, although the constitutive model was derived for dilute fiber suspensions. In the computations, a fiber was assumed to align along the streamline (co-linear fiber alignment). Lipscomb et al. (1988) stated that the continuum theory for dilute fiber suspensions that incorporated the orientation distribution function into the stress equation predicted the flow field accurately. Furthermore, Baloch and Webster (1995) conducted numerical simulations for dilute to semi-dilute fiber suspension flows in a Newtonian solvent through various contraction and expansion geometries using the constitutive model of Lipscomb. They adopted co-linear and orthogonal fiber alignment (the local fiber orientation is perpendicular to the streamline) conditions and compared the simulated results with the experimental observations of Abdul-Karem et al. (1993).

* On the other hand, the characterization of the fiber orientation in a recirculating flow is very important from a practical view point because a recirculating flow occurs in a complex flow which is often seen in polymer processing industries. Chang et al. (1994) investigated the flow-induced orientation of liquid crystalline polymers in a planar contraction, expansion and in a dimpled channel. They employed the simplified Leslie-Ericksen equation (Ericksen (1960), Leslie (1966) (1968)) with the high viscosity approximation and found the director orientation in a recirculating region. Recently Chiba et al. computed the two-dimensional fiber orientation in recirculating flows of Newtonian fluid within a slot mounted on a wall of a parallel plate channel (Chiba et al. (1996-a), (1996-b)) and within a salient corner of an abrupt expansion channel (Chiba et al. (1997)). They reported that steady and almost complete orientation was achieved after several circulations of fibers along a streamline in a recirculating flow. The fiber orientation in steady recirculating flows was also reported by Azaiez et al. (1997), Poitou et al. (2000) and Chinesta and Chaidron (2001).
In numerical simulations for fiber suspension flow and fiber orientation, the following
two types of strategy are well known:

I. The coupled equations of the governing equations for flow field, the constitutive equa-
tion and the evolution equation of the orientation tensor with a closure equation are
solved.

II. The coupled equations of the governing equations for flow field, the constitutive equa-
tion and the evolution equation of orientation distribution function, e.g., the Fokker-
Plank equation, are solved.

* Although a statistical orientation distribution is necessary to describe the ori-
entation dependent properties of fiber suspensions completely, direct solution of
this function, such as solution of the Fokker-Planck equation, or its construction
from the solution of numerous discrete fibers did not gain popularity because of
mathematical complexity and computational intensity. As a consequence, the
implementation of tensorial quantities such as orientation tensors found some
acceptance for microstructural characterization in the works of Hinch and Leal
(1976) on suspension rheology, Doi (1981) and Grmela and Carreau (1987) on
conformation tensor rheological model. Orientation tensors provided consider-
able advantages for numerical calculations because of their concise description
of the orientation state. The major drawback of the orientation tensors was the
requirement of a closure equation in order to approximate the higher order ten-
sors which appeared in the orientation evolution equations. Advani and Tucker
(1987) (1990), Altan et al. (1989) and Chiba and Nakamura (1995) tackled the
problems related to closure approximations.

* Akbar and Altan (1992) discussed the drawbacks of using the approach based on orientation tensors, and proposed a statistical method which made it possible to generate the orientation distribution function by considering large numbers of fibers without solving the Fokker-Planck equation. Hence, the solution of orientation fields along the streamlines in a complex flow was feasible by considering numerous fibers starting from specified initial orientations. This statistical method was recently used by Ericsson *et al.* (1997) in computations of the rheology of a discontinuous fiber filled polypropylene in a squeezing flow between parallel plates. It was shown that the statistical method gave an excellent prediction of the behavior of the exact solution in simple flows using only 100 test fibers. Furthermore, Chiba and Nakamura (1995) used the statistical method to examine the accuracy of the quadratic closure approximation of the fourth-order orientation tensor.

2 MECHANICAL MODELING: CONSTITUTIVE EQUATION

Two main steps allow to derive macroscopic constitutive equations for anisotropic suspen-
sions: (i) a standard homogenization procedure that leads to volume averaged quantities if
the orientation of the particles is perfectly defined, and (ii) a statistical description of the
orientation.

**Volume averaged quantities.** In this procedure, the basic equations are obtained with a volume averaged procedure (Batchelor (1970-a)). Let $\Omega$ be a representative volume of the macroscopic scale, which contains many particles located in $\Omega_i$. Let $\tau$, $\vec{u}$ and $\mathbf{D}$ denote respectively the microscopic stress tensor, velocity vector and strain rate tensor and $\vec{n}$ be a normal vector to the particle boundary (Figure 9).
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Representative volume at a macroscopic scale</figcaption></figure>

For a given fiber orientation $\vec{p}$ (unit vector aligned along the principal axis of the fiber) in the suspension, for a Newtonian ambient fluid of viscosity $\eta$ and if inertia terms can be neglected, the corresponding macroscopic variables $\mathbf{T}$, $\vec{v}$, and $\mathbf{D}$ are related by

$$ \mathbf{T}(\vec{p}) = -p\mathbf{I} + 2\mu\mathbf{D}(\vec{v}) + \sum_p \qquad (2.1) $$

$$ \sum_p = \frac{1}{\Omega} \sum_i \oint_{\partial\Omega_i} \tau \vec{n} \otimes \vec{x} \, dx \qquad (2.2) $$

In the above expressions, $p$ is the volume averaged pressure, $\mathbf{I}$ is the unit tensor and $\sum_p$ represents the contribution of the particle to the averaged extra stress tensor.

**Statistical average.** The fiber orientation $\vec{p}$ is given in an average sense by its probability distribution $\psi(\vec{p})$ so that the macroscopic statistical averaged stress tensor (i.e. the stress tensor used in engineering computations) is

$$ \boldsymbol{\sigma} = \int \mathbf{T}(\vec{p}) \psi(\vec{p}) \, d\vec{p} \qquad (2.3) $$

For most engineering applications, the exact calculation of $\int \psi(\vec{p})$ at each point is neither necessary nor possible. This orientations information is correctly approximated with the use of the second order and fourth order orientation tensors, defined as follows:

$$ \mathbf{a} = \int \vec{p} \otimes \vec{p} \psi(\vec{p}) \, d\vec{p} \qquad (2.4) \\ \mathbf{a}_4 = \int \vec{p} \otimes \vec{p} \otimes \vec{p} \otimes \vec{p} \psi(\vec{p}) \, d\vec{p} $$

where $\otimes$ denotes the tensorial product.

The flow around a single particle in a Newtonian fluid gives a first approximation for the evolution of one particle orientation. From Jeffery's calculations (Jeffery (1922)), if $\lambda$ denotes the aspect ratio of the particle, and if $\Omega$ and $\mathbf{D}$ are respectively the volume averaged vorticity and strain rate tensors, then

$$ \frac{d\vec{p}}{dt} = \Omega\vec{p} + k\{\mathbf{D}\vec{p} - (\vec{p}^T \mathbf{D}\vec{p})\vec{p}\} \qquad (2.5) $$
with

$$k = \frac{\lambda^2 - 1}{\lambda^2 + 1} \qquad (2.6)$$

The evolution equation of the probability distribution is given by

$$\frac{d\psi(\vec{p})}{dt} + \frac{\partial}{\partial \vec{p}} \left\{ \psi(\vec{p}) \frac{d\vec{p}}{dt} \right\} = \frac{\partial}{\partial \vec{p}} \left\{ D_r \frac{\partial \psi(\vec{p})}{\partial(\vec{p})} \right\} \qquad (2.7)$$

In the above equation, the rotary diffusion coefficient $D_r$ which vanishes at zero shear rates, accounts for the hydrodynamic interactions between fibers. The evolution equation for the second order orientation tensor then becomes

$$\frac{\delta \mathbf{a}}{\delta t} = \frac{d\mathbf{a}}{dt} - (\Omega \mathbf{a} - \mathbf{a}\Omega) - k(\mathbf{D}\mathbf{a} + \mathbf{a}\mathbf{D}) = -2k\mathbf{a}_4 : \mathbf{D} + D_r \left( \mathbf{a} - \frac{\mathbf{Id}}{3} \right) \qquad (2.8)$$

where the symbol “:” denotes the tensorial product twice contracted. This equation can be easily obtained considering the time derivative of Eq. (2.4) taking into account Eqs. (2.5) and (2.7).

This equation does not allow to determine the orientation tensor if the velocity field is given, because it involves the fourth order orientation tensor $\mathbf{a}_4$. The evolution equation for the fourth order tensor involves the 6th order tensor, and in general, the evolution equation for any order moment of $\vec{p}$ involves the next higher-order moment. Therefore, to obtain a closed set of equations, one needs to determine $\mathbf{a}_4$ as a function of $\mathbf{a}$. A variety of closure approximations have been proposed to relate these two orientation tensors. Among these, we will mention the quadratic closure approximation in which $\mathbf{a}_4$ is defined as:

$$\mathbf{a}_4 = \mathbf{a} \otimes \mathbf{a} \qquad (2.9)$$

## 2.1 Explicite Calculations for $\Sigma_p$

The calculation of $\Sigma_p$ as a function of the macroscopic velocity field $\vec{v}$ requires to solve a microscopic problem, which can take different forms according to the level of approximation. To date, two major approaches have been used to determine $\Sigma_p$. In the first approach based on slender bodies theory, the stress distribution at the particles boundary is reduced to a multipolar distribution. The Batchelor's model (Batchelor (1970-a)) is deduced, which has been modified by Shaqfeh and Fredrickson (1990). All of these models can be written in the following form

$$\Sigma_p = 2\eta \mathrm{Tr}[\mathbf{D}(\vec{v})(\vec{p} \otimes \vec{p})]\vec{p} \otimes \vec{p} \qquad (2.10)$$

In this expression, $\eta$ is calculated as a function of the particles concentration and of its aspect ratio which is assumed to be very large. The exact expression depends on the level of approximations. The interactions between fibers are accounted for in a certain sense, but for long or very long fibers only. A second class of models deals with ellipsoidal or spheroidal particles. This case has been extensively studied for linear elastic materials (Gilormini and Vernusse (1992), Mura (1982)). However, it turns out that it has not been so much studied for suspensions, probably because Eshelby's works (Eshelby (1957)) are not popular in the fluid mechanics scientific community.

**Suspension of spheres.** The case of a suspension of non Brownian spheres is meaningful because the suspension remains Newtonian (viscosity $\eta_{eq}$). In this case, if $\varphi$ denotes the volume fraction of spheres, and if there are no hydrodynamic interactions (dilute sus-
pension), the dilute approximation leads to Einstein's formula

$$ \Sigma_p = 5\eta\varphi\mathbf{D}(\vec{v}) \quad \Rightarrow \quad \eta_{eq} = \eta(1 + 2.5\varphi) \qquad (2.11) $$

To account for hydrodynamic interactions, a natural way consists of calculating $\Sigma_p$ by assuming that, instead of being immersed in the ambient fluid of viscosity $\eta$, the consequences of hydrodynamic interactions are summarized by placing the sphere in an ambient fluid of viscosity $\eta_{eq}$ (self-consistent model)

$$ \Sigma_p = 5\eta_{eq}\varphi\mathbf{D}(\vec{v}) \quad \Rightarrow \quad \eta_{eq} = \frac{\eta}{1 - 2.5\varphi} \qquad (2.12) $$

For small but non infinitesimal volume fractions, the self-consistent model allows to derive (with a very different approach than the original one), the Batchelor's approximation

$$ \eta_{eq} = \eta(1 - 2.5\varphi + 6.25\varphi^2) \qquad (2.13) $$

However the self-consistent model suffers from strong limitations because the equivalent viscosity increases to infinity when $\varphi = 40\%$. For this reasons, another approach is needed to improve Eq. (2.12). Following Christensen (1990), we assume that the spheres are added progressively, step by step (differential model). At each step, the suspension is macroscopically a viscous fluid of viscosity $\eta_{eq}(\varphi)$. We then add a quantity $d\varphi$ of spheres. This quantity is small enough such that Einstein's formula for dilute suspension remains valid. The only difficulty consists in noting that the reference ambient fluid is the suspension of volume fraction $\varphi$, so that the volumetric rate induced by $d\varphi$ is $d\varphi/(1 - \varphi)$

$$ d\Sigma_p = 5\eta_{eq}\frac{d\varphi}{1-\varphi}\mathbf{D}(\vec{v}) \quad \Rightarrow \quad \eta_{eq} = \frac{\eta}{(1-\varphi)^{2.5}} \qquad (2.14) $$

This differential model leads to an expression similar to that of Krieger and Dougherty (1959). However, this expression is not valid when direct contact between particles are to be considered because the global behavior of the suspension is then neither linear, nor homogeneous. In the following, we extend Eqs. (2.12) and (2.14) to ellipsoids. Two difficulties must then be overcome. The first is that a suspension of non spherical particles is non Newtonian. The second is that the calculation is more technical.

## Flow around an ellipsoidal particle immersed in an anisotropic suspension

### Algebraic preliminary result

A linear operator on the subspace of deviatoric (traceless) tensors, which exhibit a transverse isotropic symmetry of axis $\vec{p}$ depends on three parameters. A basis of such operators can be defined as $\mathbf{A}_1, \mathbf{A}_2, \mathbf{A}_3$:

$$ \left\{ \begin{array}{l} \mathbf{A}_1 = \mathrm{Id} - \mathbf{A}_2 - \mathbf{A}_3 \\ \\ \mathbf{A}_2 = \frac{1}{2} (3\vec{p} \otimes \vec{p} - \mathrm{Id}) \otimes (\vec{p} \otimes \vec{p}) \\ \\ \mathbf{A}_3 \mathbf{Q} = \mathbf{Q}(\vec{p} \otimes \vec{p}) + (\vec{p} \otimes \vec{p}) \mathbf{Q} - 2\mathrm{Tr}[\mathbf{Q}(\vec{p} \otimes \vec{p})](\vec{p} \otimes \vec{p}) \end{array} \right. \qquad (2.15) $$

This basis has two interesting properties. It exhibits the direction of transverse isotropy $\vec{p}$ and it allows easy products because it can be checked that:

$$ \mathbf{A}_i : \mathbf{A}_i = \mathbf{A}_i \quad \text{and} \quad \mathbf{A}_i : \mathbf{A}_j = 0 \quad \text{if } i \neq j \qquad (2.16) $$
## Jeffery's problem in an anisotropic fluid

Jeffery derived the exact solution for the flow of a Newtonian fluid around a particle subjected to homogeneous conditions at infinity. Eshelby (1957) gave an alternative technique, which has been generalized to the case of an anisotropic surrounding matrix. His results can be summarized as follows.

If $D(\vec{v}) = M\Sigma$ is the constitutive relation of the fluid: $M\Sigma = \frac{1}{2\eta} \Sigma$ in the isotropic case, or more generally, $M = \frac{1}{\alpha_1}\mathbf{A}_1 + \frac{1}{\alpha_2}\mathbf{A}_2 + \frac{1}{\alpha_3}\mathbf{A}_3$ in the transverse isotropic case, then, the solution to generalized Jeffery's problem reads:

$$S_{esh} M \Sigma_p = \varphi \mathbf{d}(v) \qquad (2.17)$$

where $S_{esh} = S_{esh}(\lambda, \alpha_1, \alpha_2, \alpha_3)$ is the fourth order analytical Eshelby tensor which depends on the shape of the particle and on the fluid properties only.

## Approximations

**Dilute case:** In the dilute case, $M\Sigma_p$ is calculated by neglecting the particles interaction. Eq. (2.17) is thus solved with $\alpha_1 = \alpha_2 = \alpha_3 = 2\eta$. This approximation leads to the classical dilute approximation for spheroids and generalizes Einstein relation to spheroids (Batcherlor (1970)).

**Self consistent approach:** Jeffery's problem is solved for an ambient fluid with transverse isotropic but unknown behavior, characterized by parameters $\alpha_1, \alpha_2, \alpha_3$. With the same argument as the one that yielded to Eq. (2.12), it is possible to derive the following implicit equations:

$$\left\{
\begin{aligned}
\alpha_1 &= 2\eta + \frac{\alpha_1\varphi}{2S_{1212}} \\
\alpha_2 &= 2\eta + \frac{\alpha_2\varphi}{S_{3333} - S_{3311}} \\
\alpha_3 &= 2\eta + \frac{\alpha_3\varphi}{2S_{1313}}
\end{aligned}
\right.
\qquad (2.18)$$

This set of equations is non linear but can be solved with a standard Newton Raphson algorithm.

**Differential approach:** Particles are added progressively, so that at each time, one has to deal with a dilute suspension of spheroids. However, the problem is here specific, because the suspending fluid is a suspension, which is non Newtonian. The system to solve is similar to equation (2.14):

$$\left\{
\begin{aligned}
d\alpha_1 &= \frac{\alpha_1 d\varphi}{(1 - \varphi)2S_{1212}} \\
d\alpha_2 &= \frac{\alpha_2 d\varphi}{(1 - \varphi)(S_{3333} - S_{3311})} \\
d\alpha_3 &= \frac{\alpha_3 d\varphi}{(1 - \varphi)2S_{1313}}
\end{aligned}
\right.
\qquad (2.19)$$

with initial conditions: $\alpha_i = 2\eta$ for $\varphi = 0$.
## 2.2 Results and Discussion

The numerical results have been obtained by computations (Meslin (1997)). In order to compare these results with other existing theories, the stress tensor is written in a similar form as the one introduced by Tucker (1991).

$$ \mathbf{T} = -p\mathbf{I} + 2\eta_1(\mathbf{D} + N_s((\vec{p} \otimes \vec{p})\mathbf{D} + \mathbf{D}(\vec{p} \otimes \vec{p}))) + N_p \mathrm{Tr}((\vec{p} \otimes \vec{p})\mathbf{D}(\vec{p} \otimes \vec{p})) \quad (2.20) $$

with these notations

$$ \left\{ \begin{array}{l} N_s = \frac{\alpha_3 - \alpha_1}{\alpha_1} \\ \\ 2\eta = \alpha_1 \\ \\ N_p = \frac{3(\alpha_2\alpha_1)}{2\alpha_1} - \frac{2(\alpha_3\alpha_1)}{\alpha_1} \end{array} \right. \qquad (2.21) $$

$\eta_1$ contains all the isotropic contributions to the viscosity (from both the solvent and the particles), while anisotropic contributions of the particles are represented by $N_p$ and $N_s$. However, like in the Batchelor's dilute model, we find that for both differential and self consistent schemes $N_p$ is always greater than $N_s$, even for particles of small aspect ratio. Therefore, the anisotropic contributions by the particles are essentially described by $N_p$. The product $\eta_1 N_p$ is equal to the parameter $\eta$ of the Shaqfeh's model and we compare spheroidal model to slender bodies model. Figure 10(a) shows the differences between the dilute spheroidal approximation (Batchelor spheroid), the Dinh and Armstrong's model (Dinh and Armstrong (1984)), Shaqfeh and Fredrickson's slender bodies model (Shaqfeh and Fredrickson (1990)), and our self-consistent and differential model, for three aspect ratios of the particles. It can be seen that even for an aspect ratio $\lambda$ of 350, neither Shaqfeh and Fredrickson's nor Dinh and Armstrong's model is correct because, for very low values of $\lambda$, the curves are not tangent to Batchelor's dilute approximation. This point has been already mentioned by Ranganathan and Advani (1993), who proposed a correction which accounts for the finite aspect ratio of the particles. Figures 10(b) compare the same models for various fibers shape factors. It is to be noted that as for hard spheres models (Christensen (1990)), the self-consistent approach overestimates the rheological parameters (this self-consistent model evidences also a vertical asymptote for a volumic rate less than one). The differential model and the modified Shaqfeh and Fredrickson's model for large aspect ratio lead to similar results in the semi-dilute range. Thus, in this sense the differential model generalizes Shaqfeh and Fredrickson's model in accounting for an averaged interaction between finite aspect ratio's fibers. In particular we can apply the differential model for prolate spheroids.

# 3 INJECTION MODELING

## 3.1 Mechanical Modeling

The flow model, which must be imposed in the domain $\Omega_t^f$ occupied by the fluid at time $t$, is given by the following equations:

* Equilibrium equation, which neglecting the inertia and mass terms, results in

$$ \mathrm{Div} \, \boldsymbol{\sigma} = \vec{0} \qquad (3.1) $$

where $\boldsymbol{\sigma}$ is the stress tensor.
<figure><img src="image_11.png" /><figcaption>Figure 10. Models comparison</figcaption></figure>

* Incompressibility condition

$$Div \vec{v} = 0 \qquad (3.2)$$

where $\vec{v}$ is the velocity field.

* The constitutive equation for fibers with a large aspect ratio, and where the Brownian effects are neglected, is

$$\sigma = -p\mathbf{I} + 2\eta(\mathbf{D} + N_p\mathbf{a}_4 : \mathbf{D}) \qquad (3.3)$$

where $p$ is the pressure field, $\mathbf{I}$ the unit tensor, $\mathbf{D}$ the strain rate tensor, $\mathbf{a}_4$ the fourth order orientation tensor, $\eta$ the suspension equivalent viscosity, $N_p$ a scalar parameter which depends on both the fiber concentration and the fiber aspect ratio and the symbol “:” denotes the tensorial product two times contracted, i.e. $(\mathbf{a}_4 : \mathbf{D})_{ij} = \mathbf{a}_{4ijkl}\mathbf{D}_{kl}$ (sum with respect to repeated indices).
* Definition of the orientation tensors: (already defined in Eq. (2.4))

The second order orientation tensor is defined by

$$ \mathbf{a} = \oint_{\vec{p}} \vec{p} \otimes \vec{p} \psi(\vec{p}) d\vec{p} \qquad (3.4) $$

where $\vec{p}$ is the unit vector aligned on the fiber axis direction, $\psi(\vec{p})$ the orientation probability function and the symbol $\otimes$ denotes the tensorial product, i.e. $(\vec{p} \otimes \vec{p})_{ij} = (\vec{p})_i(\vec{p})_j$.

The fourth order orientation tensor is defined in a similar way by

$$ \mathbf{a}_4 = \oint_{\vec{p}} \vec{p} \otimes \vec{p} \otimes \vec{p} \otimes \vec{p} \psi(\vec{p}) d\vec{p} \qquad (3.5) $$

* Conservation balance of the fiber orientation distribution

If we denote by $\vec{p}_0$ the orientation of a fiber located at time $t = 0$ in the point $\vec{x}_0$, then the orientation of this fiber at time $t$, when it is located in the point $\vec{x} = \vec{x}(\vec{x}_0, t)$, will be referred to as $\vec{p}$. The fiber orientation distribution at time $t$ can be deduced from the initial one from the following conservation balance

$$ |\vec{dp}_0|\psi_0(\vec{p}_0) = |\vec{dp}|\psi(\vec{p}) \qquad (3.6) $$

where $\psi_0(\vec{p})$ denotes the initial fiber orientation distribution and $|\vec{dp}_0|$ an infinitesimal solid angle defined around $\vec{p}_0$. To clarify the notation we consider a planar fiber distribution, where the fiber orientation can be defined by the angle $\varphi$ with respect to a reference axis. Thus, the previous balance can be rewritten as

$$ d\varphi_0\psi_0(\varphi_0) = d\varphi\psi(\varphi) \qquad (3.7) $$

where $d\varphi$ denotes the infinitesimal angle defined around the angle $\varphi$

$$ d\varphi = \left(\varphi + \frac{d\varphi}{2}\right) - \left(\varphi - \frac{d\varphi}{2}\right) \qquad (3.8) $$

In general the initial fiber orientation distribution is known, and as we are going to show, the final orientation $\vec{p}$ can be computed from the initial one $\vec{p}_0$, i.e. we can compute the orientation angle $\varphi$ of a fiber located at time $t$ in the point $\vec{x}$, from its initial orientation angle $\varphi_0$. Thus, we can also compute the value of $d\varphi$ from $d\varphi_0$.

The orientation probability associated with the angle $\varphi$ can be finally computed from

$$ \psi(\varphi) = \frac{d\varphi_0\psi_0(\varphi_0)}{d\varphi} \qquad (3.9) $$

* Fiber orientation transport

In order to define the relation between the final and initial orientations $\varphi = \varphi(\varphi_0, t)$ we start by introducing Jeffery's equation which governs the orientation evolution of an ellipsoidal fiber immersed in a Newtonian fluid

$$ \frac{d\vec{p}}{dt} = \Omega\vec{p} + k\mathbf{D}\vec{p} - k(\vec{p}^T\mathbf{D}\vec{p})\vec{p} \qquad (3.10) $$
where $k$ is a scalar parameter depending on the fiber aspect ratio $\lambda$ (ratio between the fiber length and diameter)

$$k = \frac{\lambda^2 - 1}{\lambda^2 + 1} \qquad (3.11)$$

and $\Omega$ represents the vorticity tensor.

The solution of the Jeffery's equation can be written by using the transport tensor $E_k$, whose evolution is governed by the following linear advection equation

$$\frac{dE_k}{dt} = (\Omega + kD)E_k \qquad (3.12)$$

and satisfies the following initial condition $E_k(t=0) = I$.

With the transport $E_k$ calculated by solving the previous advection equation, the Jeffery's solution is given by

$$\vec{p} = \frac{\mathbf{E}_k \vec{p}_0}{\|\mathbf{E}_k \vec{p}_0\|} \qquad (3.13)$$

or

$$\vec{p}_0 = \frac{\mathbf{E}_k^{-1} \vec{p}}{\|\mathbf{E}_k^{-1} \vec{p}\|} \qquad (3.14)$$

The previous equations define the flow model of a suspension of short fibers. The rheological parameters $\eta$ and $N_p$ must be identified experimentally for the high fiber concentration usually employed in industrial applications (only in the case of dilute suspensions explicit expressions for both parameters exist). In the previous model we have neglected the Brownian effects, but as discussed in other works (Meslin (1997)) the fact of considering the fibers with an ellipsoidal shape, i.e. $k \neq 1$ (even for $k$ very close to the unity) induces a continuous rotation of the fiber in the flow (even in steady shear flows). This produces an extra-dissipation that in slender body theories are induced by the introduction of an orientation diffusion term (Brownian effects). This diffusion term prevents the fully alignment of the fiber with the flow, generating the extra-dissipation experimentally observed.

**Remark:** The flow model previously established can be also applied in the simulation of steady flows, as encountered in extrusion processes. In this case special attention must be made in the treatment of steady recirculating areas where neither the initial orientation nor boundary conditions are defined. In this case only the periodicity of the solution can be imposed to find the steady solution. More details concerning the numerical strategies to compute steady solutions in steady recirculating flows can be found in Poitou et al. (2000), Chinesta and Chaidron (2001), Chaidron and Chinesta (2002-a), (2002-b), Chiba and Nakamura (1998) and Azaiez et al. (1997).

## Boundary conditions

We denote by $\Omega_t^f$ the domain occupied by the fluid at time $t$. The boundary of the fluid domain which will be referred to as $\partial\Omega_t^f$ can be divided in two parts. In the first part, denoted by $\partial_1\Omega_t^f$, the velocity is assumed known, whereas in the second one, $\partial_2\Omega_t^f$, the traction is imposed. Moreover, the boundary $\partial_1\Omega_t^f$ can be divided in other two parts: the inflow boundary $\Gamma^-$ (injection gate, which is considered fixed in time) and the wall boundaries $\Gamma_t^0$ where the adherence condition impose a zero velocity. In injection simulations,
the boundary $\partial_2\Omega_t^f$ corresponds usually to the flow front where no tractions are imposed by the empty part of the mold. The inflow boundary conditions are given by the injection flow rate and by the initial fiber orientation distribution. We can summarize the boundary conditions to impose in the flow domain at time $t$ by:

$$ \vec{v}(\vec{x} \in \partial_1 \Omega_t^f) = \begin{cases} \vec{v}_i & \vec{x} \in \Gamma^- \\ \vec{0} & \vec{x} \in \Gamma_t^0 \end{cases} \qquad (3.15) $$

$$ \sigma \vec{n}(\vec{x} \in \partial_2 \Omega_t^f) = \vec{0} \qquad (3.16) $$

$$ \psi_0(\vec{p}_0)(\vec{x} \in \Gamma^-) = \psi_i(\vec{p}_0) \qquad (3.17) $$

## 3.2 Fluid Domain Evolution: Flow Front Tracking

As commented, the flow model is defined in the part of the whole domain (injection mold) $\Omega$ occupied by the fluid at each time, $t$, $\Omega_t^f$. In order to update the fluid domain we are going to introduce the fluid presence function, also known as fluid fraction function $I(\vec{x}, t)$. This function takes a unit value in the fluid region and it is zero in the empty domain $\Omega_t^e = \Omega - \Omega_t^f$:

$$ I(\vec{x}, \in \Omega, t) = \begin{cases} 1 & \vec{x} \in \Omega_t^f \\ 0 & \vec{x} \in \Omega_t^e \end{cases} \qquad (3.18) $$

The evolution of this function is given by the following scalar and linear advection equation

$$ \frac{dI}{dt} = \frac{\partial I}{\partial t} + \vec{v} \text{Grad} I = 0 \qquad (3.19) $$

This equation is defined in the whole domain. The fluid presence function must verify a boundary condition on the mold inlet (inflow boundary)

$$ I(\vec{x} \in \Gamma^-, t) = 1 \qquad (3.20) $$

as well as an initial condition. If we assume the mold empty at the beginning of the injection process the initial condition is given by

$$ I(\vec{x} \in \Omega, t = 0) = 0 \qquad (3.21) $$

## 3.3 Numerical Modeling of the Filling Process

In the simulation of the injection process we consider a fully Eulerian description combined with an explicit strategy, solving the flow kinematics for a given fluid domain, in which the fiber orientation is assumed known. With the flow kinematics just computed we can update the fluid domain as well as the fiber orientation within the new fluid domain.

This explicit strategy allows to get an accurate treatment of the equations of motion defining an elliptic anisotropic Stokes problem (Eqs. (3.1), (3.2) and (3.3)) using a standard velocity-pressure mixed finite element formulation, as well as of the advection equations governing the fluid domain and fiber orientation evolutions. The hyperbolic character of advection equations requires specific numerical techniques for its discretisation. In this paper, the discontinuous finite element method (also known as Lesaint-Raviart technique) will
be applied to solve advection equations (Eqs. (3.12) and (3.19)). Moreover the integration of both advection problems may be carried out simultaneously.

The simulation process finishes when the mould domain is full filled by the fluid, i.e. when $\Omega_t^f = \Omega$.

### 3.3.1 Discretisation technique: equations of motion

The discretisation of the equations of motion is carried out by means of a standard mixed finite element technique, using a $P_1$ + “bubble” in the velocity interpolation and $P_1$ in the pressure approximation. This functional approximation verifies the LBB condition (Pironneau (1989)). In order to extend the formulation to the whole mold domain, we impose a pseudo-behavior in the empty volume, defined by $p = 0$ and $\vec{v} = \vec{0}$ (Pichelin and Coupez (1998)). Multiplying both variational formulations by a function of the fluid fraction $f(I)$ and $1 - f(I)$ respectively we obtain the variational formulation extended to the whole domain:

Find $\vec{v} \in (H^1(\Omega))^3$ and $p \in L^2(\Omega)$ verifying the boundary conditions on $\partial\Omega: \vec{v}(\vec{x} \in \Gamma^-) = \vec{v}_i$ and $\vec{v}(\vec{x} \in \partial\Omega - \Gamma^-) = \vec{0}$, such that

$$ 
\begin{aligned}
& \int_{\Omega} f(I) \boldsymbol{\sigma} : \mathbf{D}^* d\Omega + \int_{\Omega} \alpha(I)(1 - f(I))\vec{v}\vec{v}^* d\Omega = 0 \\
& \int_{\Omega} f(I) \mathrm{Div} \, \vec{v} p^* d\Omega + \int_{\Omega} \beta(I)(1 - f(I))p p^* d\Omega = 0
\end{aligned}
\qquad (3.22) 
$$

$\forall \vec{v}^* \in (H_0^1(\Omega))^3$ and $\forall p^* \in L^2(\Omega)$, and where the expression of the stress tensor is given by the constitutive equation (3.3).

In the previous variational formulation $H^1(\Omega)$ and $L^2(\Omega)$ denote the usual Sobolev and Lebesgue functional spaces and $H_0^1(\Omega)$ is the functional space of velocities in $H^1(\Omega)$ vanishing on the domain boundary $\partial\Omega$. The election of the functions $f(I)$, $\alpha(I)$ and $\beta(I)$ is a key point to obtain a numerical scheme without numerical dissipation and with a low diffusion of the flow front. The necessity to introduce functions $\alpha(I)$ and $\beta(I)$ follows from the dimensional consistency requirements of both integrals in each equation.

Pichelin and Coupez (1998) considered a linear combination of both variational formulations, i.e. the fluid fraction was taken as weight function, i.e. $f(I) = I$. Both constant parameters $\alpha$ and $\beta$ were considered in that work. To justify the required values of both parameters, we consider an element $\Omega^e$ which just starts its filling process. In this case the fluid fraction in this element takes an intermediate value $0 < I^e \ll 1$. If the parameters $\alpha$ and $\beta$ are too high, the over-imposition of the zero velocity condition in the empty region or in the elements which start their filling process, tends to make the natural flow front movement difficult. On the other hand, if the parameters $\alpha$ and $\beta$ are too small, the equation of motion as well as the fluid incompressibility will be over-imposed in the element $\Omega^e$ in spite of the little amount of fluid existing in this element which just starts its filling process. In this way, this pseudo-incompressibility, derived from the small values of both parameters, retards the filling process. In the same way, the resulting flow front thickness increases, and the localization of the flow front remains uncertain. In order to improve the flow front localization (reducing its numerical thickness, i.e. the number of elements in the flow direction with a fluid fraction $0 < I < 1$) Coupez et al. (1999) have proposed to use small values of both parameters combined with a local mesh adaptation without topological changes. This mesh adaptation tends to reduce the volume of the elements located on the flow front ($0 < I < 1$).
Another method able to reduce the flow front thickness was proposed by Monton (2002). In this case the function $f(I)$ proposed is given by

$$f(I^e) = \begin{cases} 1 & I^e \ge I_{th} \\ 0 & I^e < I_{th} \end{cases} \qquad (3.23)$$

where $I_{th}$ represent a threshold value, close to one. The first advantage of this election is that no combinations of the fluid flow and the pseudo-behavior variational formulations are employed in the partially filled elements. Thus, only in the case that $I^e \ge I_{th}$ the flow model will be enforced in the element $\Omega^e$. However, as we shall show later, a judicious choice of the parameters $\alpha$ and $\beta$ will is essential to enforce the mass conservation reducing the flow front thickness to one element.

To illustrate the parameters choice, we consider a one-dimensional domain. We denote by $\Omega^e$ a generic element where a linear approximation of the velocity field is considered. Thus, the velocity in each point of $\Omega^e$ can be expressed from its nodal values $v_e$ and $v_{e+1}$. The coordinates of both nodes are $x_e$ and $x_{e+1}$ ($x_{e+1} > x_e \ge 0$) respectively. We assume the injection gate located at $x = 0$, where the injection velocity is enforced to a value $v_i$. Now, let $I^e(t) < I_{th}$ be the fluid fraction existing in the element $\Omega^e$ at time $t$. We will assume that the element located upstream from $\Omega^e$, $\Omega^{e-1}$, is fully filled, i.e. $I^{e-1}(t) = 1$ and the downstream element remains empty, i.e. $I^{e+1}(t) = 0$.

Thus, at time $t$, with $I^{e-1}(t) = 1$, the flow equations are imposed in the element $\Omega^{e-1}$. The velocity at node $x_e$ tends to the injection velocity as a direct consequence of the fluid incompressibility and the one-dimensional flow. On the other hand, the node $x_e$ also belongs to the element $\Omega^e$ where due to the fact that $I^e(t) < I_{th}$ a null velocity is imposed, which implies the annulations of both nodal velocities, and in particular the annulation of the velocity at node $x_e$. An evident conflict appears because each element containing the node $x_e$ tends to impose a different value of the velocity at this node. Due to the continuity assumed in the velocity interpolation the resulting velocity in this node will be an intermediate value between the injection velocity (required by the flow equations) and a null velocity (required by the pseudo-behavior). Therefore, the extended variational formulation to the whole domain, as has been presented until now, affects the natural flow front movement. To avoid these undesired effects one possibility consists in a clever choice of the parameters $\alpha$ and $\beta$ (Monton (2002)). The main idea is to enforce the imposition of the flow equations with respect to the pseudo-behavior model. If we consider the following parameters

$$\alpha(I^e) = \begin{cases} 1 & I^e = 0 \\ \varepsilon_\alpha & I^e > 0 \end{cases} \qquad (3.24)$$

and

$$\beta(I^e) = \begin{cases} 1 & I^e = 0 \\ \varepsilon_\beta & I^e > 0 \end{cases} \qquad (3.25)$$

where $\varepsilon_\alpha$ and $\varepsilon_\beta$ are two constants small enough, the nodal velocity at node $x_e$ is dominated by the flow equations imposed in the element $\Omega^{e-1}$ (the contribution of the pseudo-behavior imposed in the element $\Omega^e(0 < I^e < I_{th})$, and consequently affected by the constants $\varepsilon_\alpha$ and $\varepsilon_\beta$ is negligible). Thus the velocity at node $x_e$ approaches the exact value (the injection velocity). The velocity at node $x_{e+1}$ will be null due to the imposition of the pseudo-behavior (null velocity) in all the elements containing this node.

**Remark.** A similar and alternative possibility consists in assembling only the contributions of elements whose fluid fractions are higher than the threshold value, i.e. $I^e > I_{th}$. The main problem of this technique is the fact that in the global equations system, the rows
related to the nodes which are connected only to empty elements remain null. In order to avoid the matrix singularity we can proceed to modify the system adding, for example, a unit value in the diagonal component of each null row. This solution is equivalent to imposing zero velocities and pressures in the nodes within the empty region. We notice that this technique is very close to the standard control volume technique widely used in the simulation of injection processes.

The proposed strategy reduces the flow front thickness to one element with a great accuracy in the mass conservation. However, the exact flow front position remains uncertain, which makes difficult, for example, the imposition of surface tension. In the simulation of usual polymers or molten composites forming processes (injection, extrusion, etc.) the surface tension can be neglected and the proposed strategy brings accurate results. Polymer processes involving multiphase such as gas assisted injection, or multi-materials flows such as co-injection and co-extrusion, require an interface tracking. Our strategy can be generalized for the treatment of these problems by the introduction of a presence function for each material and for each phase. The evolution of these presence functions can be computed by solving the corresponding advection equations governing their evolutions. However, as commented previously, the interface location remains lightly uncertain. Another strategy able to treat accurately interface problems is the level set method (Sussman et al. (1994)). This technique uses a level function $\phi$ whose zero value determines the interface position. The evolution of this function is given by the following linear advection equation

$$ \frac{\partial \phi}{\partial t} + \vec{u} \text{Grad} \phi = 0 \qquad (3.26) $$

where $\vec{u}$ is an arbitrary velocity field which must coincide with the material velocity $\vec{v}$ on the interface at time $t$, and which will denoted from now on by $\Gamma_i(t)$:

$$ \vec{u}(\vec{x} \in \Gamma_i(t), t) = \vec{v}(\vec{x} \in \Gamma_i(t), t) \qquad (3.27) $$

It is easy to verify that the solution of Eq. (3.26) allows to update the interface. Thus, if we assume that the level function at time $t$, $\phi(\vec{x}, t)$, as well as the flow kinematics $\vec{v}(\vec{x}, t)$ are known, then the interface at this time is given by

$$ \Gamma_i(t) = \{\vec{x}/\phi(\vec{x}, t) = 0\} \qquad (3.28) $$

The solution of Eq. (3.26) allows to compute the level function at time $t + \Delta t$, $\phi(\vec{x}, t + \Delta t)$, from which we can update the interface position

$$ \Gamma_i(t + \Delta t) = \{\vec{x}/\phi(\vec{x}, t + \Delta t) = 0\} \qquad (3.29) $$

In the above equation, $\Delta t$ is the time step used in the time discretisation of Eq. (3.26) that we will describe in detail more latter. The main difficulty of this technique is that the quality of function $\phi(\vec{x}, t)$ is degraded during the process simulation. Thus, the slope of the level function can develop into very high or very small values. In both cases a significant error can appear in the interface updating described in the previous paragraphs. In order to reduce these numerical errors some researchers proposed to modify at each iteration the level function in order to obtain a unit slope in the flow direction on the interface. To this purpose, and knowing the level function at time $t$, we need to solve the following non linear advection equation

$$ \frac{\partial \bar{\phi}}{\partial t} + \operatorname{sgn}(\phi) \frac{\text{Grad} \bar{\phi}}{\|\text{Grad} \bar{\phi}\|} \text{Grad} \bar{\phi} = \operatorname{sgn}(\phi) \qquad (3.30) $$
where the function $\bar{\phi}$ coincides with $\phi$ along the interface, i.e.

$$ \bar{\phi}(\vec{x} \in \Gamma_i(t), t) = \phi(\vec{x} \in \Gamma_i(t), t) \qquad (3.31) $$

and it has an unit slope on the interface in the flow direction.

The numerical treatment of non linear advection equations is a difficult matter. However the main difficulty in the use of the level set strategy is the necessity of having an accurate velocity description on the interface. To attain this accuracy, the consideration of an incompressible pseudo fluid filling the empty region of the mould is a good alternative. Thus, usually, a virtual incompressible fluid, with a very low viscosity (a gas for example), is assumed to occupy the empty region of the mould. In this case, a continuous velocity field is obtained, which allows an accurate determination of the velocity on the interface. However, the fact of using an incompressible pseudo fluid enforces the disposition of some holes on the mold walls to evacuate this pseudo fluid, allowing the flow front movement of the injected fluid. When an evacuation hole is reached by the fluid, we need to obstruct the hole to avoid fluid losses. On the other hand, if the fluid reaches a region without holes, the incompressibility of the pseudo fluid makes impossible the filling of this region. The optimal holes position, their number, opening control, etc. require some additional considerations. Another difficulty results from the low viscosity of the pseudo fluid which generates boundary layers in the velocity field in the neighboring of the walls due to the adherence (stick) boundary condition characteristic of the viscous fluids. Some researchers propose the substitution of the stick boundary condition by a slip boundary condition in the region occupied by the pseudo fluid. In this case, an adequate control in the imposition of the boundary conditions is required. Some references concerning the use of the level set techniques in the context of polymer injection forming processes can be found in VerWeyst et al. (1999). However, works concerning the efficient treatment of the numerical difficulties previously described are rare.

For all these considerations and without a precise interface location requirement, we have preferred to use the extended variational formulation coupled with a fluid presence function description.

### 3.3.2 Discretisation technique: advection equations

In this section we will describe an accurate numerical treatment of the advection equations. In order to alleviate the notation we consider the linear and scalar advection equation governing the evolution of the fluid presence function $I$

$$ \frac{dI}{dt} = \frac{\partial I}{\partial t} + \vec{v} \text{Grad } I = 0 \qquad (3.32) $$

Taking into account the fluid incompressibility, $Div \vec{v} = 0$, Eq (3.32) can be rewritten in the form

$$ \frac{dI}{dt} = \frac{\partial I}{\partial t} + Div(\vec{v}I) = 0 \qquad (3.33) $$

Now, we can consider the integral conservation form related to the previous differential equation

$$ \int_{\Omega^e} \frac{\partial I}{\partial t} d\Omega + \int_{\partial \Omega^e} Div(\vec{v}I) d\Omega = 0 \qquad (3.34) $$
for any subdomain, that we can identify with the elements $\Omega^e$ of a finite element mesh. Using the divergence theorem, the previous equation becomes

$$ \int_{\Omega^E} \frac{\partial I}{\partial t} d\Omega + \int_{\partial \Omega^e} \vec{v} \vec{n} I \, dS = 0 \qquad (3.35) $$

where $\partial \Omega^e$ denotes the boundary of $\Omega^e$, and $\vec{n}$ is the outwards unit vector defined on the element boundary. The element boundary can be divided in two parts: $\partial \Omega_+^e$ through which the flow leaves the element $\Omega^e$ (i.e. the part of the element boundary verifying the relation $\vec{v}\vec{n} > 0$); and the inflow element boundary $\partial \Omega_-^e$, though which the fluid is coming to the element $\Omega^e$ (verifying $\vec{v}\vec{n} < 0$). Thus, the mass balance results

$$ \int_{\Omega^E} \frac{\partial I}{\partial t} d\Omega + \int_{\partial \Omega_-^+} \vec{v} \vec{n} I \, dS + \int_{\partial \Omega_+^e} \vec{v} \vec{n} I \, dS = 0 \qquad (3.36) $$

Now we can introduce a constant approximation of the fluid fraction $I$ into each element. In this case the fluid fraction is not defined in the element boundary, and the flow coming in the element $\Omega^e$ is approximated taking into account the fluid fraction existing in the upstream element $\Omega^{e-}$, $I^{e-}$. In the same way, the outflow will depend on the fluid fraction existing in the element $\Omega^e$. Thus, the previous equation can be written as

$$ \frac{\partial I^e}{\partial t} |_{\Omega^e} = \sum_i q_i^{e-} I_i^{e-} - q_i^{e+} I^e \qquad (3.37) $$

where $|\Omega^e|$ is the volume of the element $\Omega^e$, $q^e = \int_{\partial \Omega_+^e} \vec{v} \vec{n} \, dS$, $q^e I^e$ represents the flow rate leaving the element $\Omega^e$, and $\sum_i q_i^{e-} I_i^{e-}$ denotes the flow coming to the element $\Omega^e$ from all the upstream elements ($I_i^{e-}$ representing the fluid fraction existing in each upstream element).

Now, we are going to treat the time derivative. The Taylor's expansion of $I^e(t + \Delta t)$ is

$$ I^e(t + \Delta t) = I^e(t) + \left. \frac{\partial I^e}{\partial t} \right|_t \Delta t + \left. \frac{1}{2} \frac{\partial^2 I^e}{\partial t^2} \right|_t (\Delta t)^2 + \dots \qquad (3.38) $$

If we consider a first order explicit scheme, the fluid fraction actualization results

$$ I^e(t + \Delta t) \approx I^e(t) + \left. \frac{\partial I^e}{\partial t} \right|_t \Delta t = I^e(t) + \frac{\sum_i q_i^{e-} - q^e I^e}{|\Omega^e|} \Delta t \qquad (3.39) $$

**Remark.** Higher order discretization schemes can be obtained directly from the Taylor's expansion which requires higher order derivatives of the fluid fraction. For their evaluation we can proceed by derivation of the mass balance

$$ \frac{\partial^p I}{\partial t^p} + \mathrm{Div} \left( \vec{v} \frac{\partial^{p-1} I}{\partial t^{p-1}} \right) = 0 \qquad (3.40) $$

Thus, operating as before we obtain

$$ \frac{\partial^p I^e}{\partial t^p} |_{\Omega^e} = \sum_i q_i^{e-} \frac{\partial^{p-1} I_i^{e-}}{\partial t^{p-1}} - q_i^{e+} \frac{\partial^{p-1} I^e}{\partial t^{p-1}} \qquad (3.41) $$
and the fluid fraction updating results

$$I^e(t + \Delta t) \approx I^e(t) + \frac{\sum_i q_i^{e-} I_i^{e-} - q^e I^e}{|\Omega^e|} \Delta t + \dots + \frac{1}{n!} \frac{\sum_i q_i^{e-} \frac{\partial^n I_i^{e-}}{\partial t^n} - q^e \frac{\partial^n I^e}{\partial t^n}}{|\Omega^e|} (\Delta t)^n \quad (3.42)$$

The treatment of the others advection equations is carried out in a similar way. The main difficulty found using this Eulerian integration strategy is related to field initializations. Thus, when an element $\Omega^e$ starts its filling process at time $t$, the fluid fraction at this time is well defined ($I^e(t) = 0$). However, if we assume a general field $T$ whose evolution is governed by a similar advection equation

$$\frac{dT}{dt} = \frac{\partial T}{\partial t} + \vec{v} \text{Grad} T = 0 \quad (3.43)$$

its value in the empty elements cannot be defined. In this case the field updating given by Eq. (3.40) remains indeterminate. A solution to this problem is proposed in Chinesta, Mabrouki and Ramon (2002), that in spite of its Eulerian character, avoid the necessity of a field initialization. In this case, when an element starts its filling process, the upstream element value governs the element updating. In this form, the proposed technique seems not longer to Lagrangian techniques, such as the method of characteristics.

## 3.4 Numerical Results and Experimental Validation

Rectangular plates of 10 × 10 × 1 cm. have been injected using PA 6-6 reinforced with short glass fibers. The viscosity and the parameter $N_p$ were identified using rheological devices and techniques proposed in Meslin et al. (1998), (1999). With these values the numerical model was solved and a prediction of the fiber orientation was obtained. Figure 11 shows a schematic of the injected pieces and two optical micrographs taken on the section shown in the figure. A good agreement between experimental and numerical results is obtained. In the graphical representation of the fiber orientation we have depicted the eigenvector associated with the higher eigenvalue, that represents the most probable fiber orientation. We can see that as expected the fibers are aligned in the flow direction in the regions of high shear rates, and in the ortho-radial direction in the neighboring of the middle plane, where a high elongation is found. Figure 12 depicts a filling sequence, as well as different views of the filling process of a 3D structural piece. We would like to point out the flow symmetry during the filling, in spite of the anisotropic viscosity, that usually gives rise to asymmetric flows. The less confined is the flow, the more relevant is this anisotropic effect (see for example Figures 12 and 14).

# 4 NUMERICAL SIMULATION AND OBSERVATION OF FIBER SUSPENSION FLOW AND FIBER ORIENTATION THROUGH A COMPLEX CHANNEL

## 4.1 Introduction

Suspensions of rodlike particles in either Newtonian or non-Newtonian fluids are often encountered in both engineering and, physical and biological sciences. The characterization of the flow or the fiber orientation in a short fiber suspension is of particular interest in current polymer processing research. When fibers form a flow-induced orientation state, fiber suspension exhibits anisotropy. Therefore, the addition of short fibers to a Newtonian liquid can drastically change the flow kinematics even at very low concentrations and in recent years the importance of coupling flow field and fiber orientation has been recognized.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Experimental validation of the fiber orientation distribution in injected plates</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Mold filling simulation of a 3D structural piece</figcaption></figure>

Furthermore, an orientation evolution model and a rheological constitutive equation that relates the stress in fiber suspensions to the orientation state are required and should be solved simultaneously to understand and describe the fiber suspension flows under certain conditions.

This section describes “fiber suspension flow and fiber orientation through a complex channel” on the basis of both numerical simulation and visualizations of flow and fiber orientation. Fiber suspensions focused in this article are rigid fiber suspensions in Newtonian liquids. The contents of this section consists of four subsections as follows:

1. Entry flow of dilute fiber suspensions through a circular contraction (Chiba *et al.* (1994-b)):

Comparison between flow predictions using co-linear alignment of fibers and flow visualizations are presented.
<figure><img src="image_14.png" /><figcaption>Figure 13. Fiber orientation evolution during the filling process of a 3D structural piece</figcaption></figure>

2. Discussion on a quadratic closure approximation using a statistical scheme (Chiba and Nakamura (1995)):
Examination of a quadratic closure approximation and the usefulness of a statistical scheme are described.

3. Analysis of fiber orientation in a recirculating flow (Chiba and Nakamura (1998)):
Computations and observations of fiber alignment in a recirculating flow within a salient corner of a 1:4 backward-facing step channel are compared.

4. Numerical solution of fiber suspension flow through a rectangular channel (Chiba et al. (2001)):
Flow and stress fields, which were computed by coupling flow kinematics with fiber orientation distribution, are discussed for the planar orientation of fibers and for the three-dimensional orientation case.
<figure><img src="image_15.png" /><figcaption>Figure 14. Schematic diagram of the fiber orientations in a fluid element flowing along a streamline</figcaption></figure>

## 4.2 Entry Flow of Dilute Fiber Suspension Through a Circular Contraction

### 4.2.1 Governing equations

The motion of an incompressible material is governed by mass and momentum conservation laws:

$$Div \vec{v} = 0 \qquad (4.1)$$

$$\rho d\vec{v}/dt = -Grad p + Div \tau \qquad (4.2)$$

where $\vec{v}$ is the velocity vector, $\rho$ is the density of fluid, $p$ is the isotropic pressure, $\tau$ is the extra stress tensor and $d/dt$ denotes the material derivative.

The continuum theory for dilute suspensions of rigid, high aspect-ratio particles developed by Lipscomb et al. (1988) is used in order to account for the presence of fibers on the stress field. The extra stress may be written as

$$\tau = 2\eta D + \phi \mu D : a_4 \qquad (4.3)$$

where $\eta$ is the viscosity of the suspending Newtonian fluid, $D$ is the deformation rate tensor, $\phi$ is the volume fraction of fibers, $\mu$ is a material constant and $a_4$ denotes the fourth-order orientation tensor, defined by Eq. (3.5).

The orientations of fibers in a fluid element, which flows along a streamline, can be distributed as illustrated in Figure 14. In the limit of high aspect-ratio of fibers, $\mu$ reduces to

$$\mu = \frac{\eta \lambda^2}{\ln(\lambda)} \qquad (4.4)$$

where $\lambda$ is the fiber aspect ratio (length to diameter ratio).

The computations of fiber suspension flows through a 4:1 circular contraction (contraction ratio $\beta = 4$) were carried out under a co-linear alignment approximation: fibers are assumed to be fully aligned with the streamlines. Therefore, the extra stress of the suspension could be simply determined from the direction of the tangent to the streamline at the location of a single fiber, i.e. the orientation unit vector $\vec{p} = \vec{v}/||\vec{v}||$ and therefore, it is not necessary to solve for the orientation distribution of fibers. The interactions of the flow field and planar fiber orientation fields were studied. For the computations, the single parameter in the extra stress equation, $\phi\mu/\eta = \phi\lambda^2/\ln(\lambda)$, was varied, which can
be determined completely by the volume fraction and aspect ratio of the fibers, $\phi$ and $\lambda$, respectively. It is called the fiber parameter.

### 4.2.2 Change of the flow patterns

The circular tube used in the flow visualizations of fiber suspension was 32 mm in diameter and 500 mm in length with a 4:1 contraction. The test fluids were seeded with polystyrene particles smaller than 0.1 mm in diameter and planar sheet of light illuminated the flow in the central plane of the tube.

The two kinds of vinylon fibers used were VF1 of diameter $d = 13~\mu\text{m}$ and length $l = 3$ mm resulting in an aspect ratio $\lambda = 231$, and VF2 of diameter $d = 11~\mu\text{m}$ and length $l = 1$ mm leading to $\lambda = 94$. These fibers were suspended in 65 wt% corn syrup/water solution of viscosity was $\eta = 0.08~\text{Pa} \cdot \text{s}$ at $16.5~^{\circ}\text{C}$.

The computed streamlines and observed flow patterns at low Reynolds number are shown in Figure 15 for a Newtonian fluid and the fiber suspension with $\phi\mu/\eta = 8$. The photograph on the right-hand side in Figure 15(b) is the flow pattern for VF1 suspension. The most obvious result is that the salient corner vortex significantly grows for the fiber suspension, despite the fact that $\phi\mu/\eta \approx 8$ corresponds to a very low volume fraction of fibers, 0.082 % for VF1 suspension. Even for large growth of the corner vortex, the vortex boundary and streamlines in the main flow become rather straight for rigid fiber suspensions. These flow patterns present striking contrast to the wine-glass shaped flow patterns with the convex vortex boundary in flexible molecule systems, such as polymer solutions (see Figure 6 in Chiba et al. (1994-b)).

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Computed streamlines and observed flow patterns through a 4:1 circular contraction at low Reynolds number; (a) Newtonian fluid; (b) fiber suspension with $\phi\mu/\eta = 8$</figcaption></figure>

Furthermore, Figure 16 shows the fiber parameter $\phi\mu/\eta$ vs. the non-dimensional vortex length $L_v^*$, where $L_v^*$ is defined as the ratio of the vortex length to the upstream tube diameter $L_v/d$. The salient corner vortex increases with an increase in $\phi\mu/\eta$, i.e. the corner vortex grows as the volume fraction and/or aspect ratio of fibers increase. This trend is in good agreement with the observations and computations. On the other hand, it can be supposed that the difference of the vortex length between VF1 and VF2 suspensions
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Comparison of vortex enhancement in a 4:1 circular entry flow of fiber suspension between computations and visualizations at $Re = 1$; computations (solid line); visualizations for VF1 suspensions (solid line with triangle); visualizations for VF2 suspensions (solid line with circle)</figcaption></figure>

is mainly due to the difference of the fiber concentration distribution in a tube: depletion region of fibers near the tube wall, especially within the salient corner vortex, becomes wider for suspensions of the longer fibers, VF1. In the computations, the fibers are assumed to be uniformly distributed, but in the experiments their distribution is not homogeneous, and the growth of the salient corner vortex strongly depends on the concentration distribution of fibers. One may also suspect that the fiber-fiber and fiber-wall interactions are significant within the corner vortex region, and can affect the flow dynamics.

The effect of the Reynolds number on the vortex length can be clearly seen in Figure 17. The vortex length is almost independent of the Reynolds number under $Re$ less than 0.1, and decreases as the Reynolds number is further increased in rigid fiber suspension flows: the main flow pushes down the salient corner vortex toward the corner and streamlines in the main flow become concave due to fluid inertia (see Figure 8 in Chiba *et al.* (1994-b)).

### 4.2.3 Stress distributions

Figure 18 shows the distribution of $T_{zz}^* - T_{rr}^*$ (extra-stress tensor components) at low Reynolds number: (a) is the distribution in Newtonian flow and (b) is the distribution in the flow of fiber suspension with $\phi\mu/\eta = 8$. Where the z-axis and r-axis are the centerline and radial directions, respectively, furthermore, the stress components are normalized such that $T_{ij}^* = T_{ij}/(2\eta U/D)$. In entry flow, the fluid element is largely stretched along the centerline as the contraction inlet is approached, and $T_{zz}^* - T_{rr}^*$ reaches its maximum at the immediate upstream region of the contraction inlet for both Newtonian and fiber suspension flows. However, the maximum value of $T_{zz}^* - T_{rr}^*$ in the Newtonian flow is about 200, while it is approximately 700 in the fiber flow. The addition of short fibers to a Newtonian liquid can drastically change the normal stress difference field even at very low concentrations. On the other hand, the shear stress field in fiber suspension flows was similar to that in
<figure><img src="image_18.png" /><figcaption>Figure 17. Computed results for the dependence of vortex length on the Reynolds number in a 4:1 circular entry flow of fiber suspension</figcaption></figure>

Newtonian case (see Figures 13 and 14 in Chiba et al. (1994)).

### 4.2.4 Conclusions

Numerical simulations and observations were presented for the flow of dilute suspensions of rigid, high aspect-ratio fibers in Newtonian fluids through a 4:1 circular contraction. The most obvious result is that, for the fiber suspension, the salient corner vortex grows significantly and the normal stress difference $T_{zz}^* - T_{rr}^*$ increases a lot even at a very low volume fraction. Furthermore, the main flow patterns become nearly conical and the corner vortex length is almost independent of flow rate for flows of rigid fiber suspensions with Reynolds numbers less than 0.1. These results present striking contrast to those for flexible molecule systems.

## 4.3 Discussion on a Quadratic Closure Approximation Using a Statistical scheme

### 4.3.1 Governing equations

The evolution equation of the orientation unit vector $\vec{p}$ for a fiber in a Newtonian flow can be described by the Jeffery equation (Eq. (3.10)) under the following assumptions: (i) the fiber may be represented by an ellipsoid of revolution; (ii) the velocity field is only locally perturbed by the motion of the fiber; (iii) the flow far from the fiber is homogeneous on a length scale that is large compared to the fiber dimensions; (iv) the motion is sufficiently slow that inertia forces are negligible; and (v) the fiber translates with the fluid velocity

$$ \frac{d\vec{p}}{dt} = \Omega\vec{p} + k(\mathbf{D} - \vec{p}^T\mathbf{D}\vec{p})\vec{p} \quad (4.5) $$

The parameter $k$ is given by

$$ k = (\lambda^2 - 1)/(\lambda^2 + 1) \quad (4.6) $$
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Distribution of the normal stress difference $T_{zz}^* - T_{rr}^*$ in a 4:1 circular entry flow at $Re = 0.01$. (a) Newtonian fluid; (b) fiber suspension with $\phi\mu/\eta = 8$</figcaption></figure>

where $\lambda$ is the fiber aspect ratio. The second-order orientation tensor $\mathbf{a}$ defined by Eq. (3.4) is a suitable and concise way of describing the orientation distribution of fibers

$$ \mathbf{a} = \oint \vec{p} \otimes \vec{p} \psi(\vec{p}) d\vec{p} \qquad (4.7) $$

then the orientation state can be illustrated using an orientation ellipse shown in Figure 19(b). The eigenvalues and the eigenvectors of $\mathbf{a}$ give the lengths and directions of the two major axes of the ellipse, and also the eigenvalues indicate the degree of orientation along these directions. $\mathbf{a}$ is symmetric and its trace is equal to unity.

On the other hand, the evolution equation of the second-order orientation tensor can be written as

$$ \frac{da}{dt} = \Omega \mathbf{a} - \mathbf{a}\Omega + \kappa(\mathbf{D}\mathbf{a} + \mathbf{a}\mathbf{D} - 2\mathbf{D} : \mathbf{a}_4) \qquad (4.8) $$

This equation involves the fourth-order orientation tensor $\mathbf{a}_4$, then a closure equation is required to approximate $\mathbf{a}_4$ using $\mathbf{a}$. We used a quadratic closure approximation given by

$$ \mathbf{a}_4 = \mathbf{a} \otimes \mathbf{a} \qquad (4.9) $$

In this section planar orientation of fibers in a developing flow of a Newtonian fluid through a parallel plate channel was analyzed by two methods to examine the performance of a quadratic closure approximation: the first method is a solution of the evolution equation of the second-order orientation tensor with a quadratic closure equation. The material
derivative $d/dt$ in Eq. (4.8) corresponds to the partial derivative $\partial/\partial t$ in a Lagrangian sense, then the evolution equations of the components of **a** may be written as

$$ \frac{\partial a_{11}}{\partial t} = 2\Omega_{12} + a_{12} + 2k(D_{11}a_{11} + D_{12}a_{12} - D_{11}a_{11}^2 - 2D_{12}a_{11}a_{12} - D_{22}a_{11}a_{22}) \quad (4.10) $$

$$ \frac{\partial a_{12}}{\partial t} = \Omega_{12}(a_{22} - a_{11}) + k(D_{12} - 2D_{11}a_{11}a_{12} - 4D_{12}a_{12}^2 - 2D_{22}a_{12}a_{22}) \quad (4.11) $$

The second method is a statistical scheme: orientation of a large number of fibers ($N = 3600$) with aspect ratio $\lambda = 5$ and 50, with specified initial angles, could be evaluated by computing the Jeffery equation along the streamlines using a 4th-order Runge-Kutta method. Then, the components of **a** could be calculated by Eq. (4.12). In this equation the orientation angle $\beta$ is an angle between the x-axis and the fiber axis which is measured counterclockwise.

$$ \begin{aligned} a_{11} &= \sum_{i=1}^{N} \cos^2 \beta_i / N \\ a_{12} &= \sum_{i=1}^{N} \sin \beta_i \cos \beta_i / N \end{aligned} \quad (4.12) $$

We examined the effect of the number of fibers used in the statistical scheme on the orientation distribution of fibers for planar orientation (Chiba and Nakamura (1994-a)) and it was found that the number of fibers of 1800 or more was appropriate. Furthermore, the orientation distribution function of fibers and the preferred angle of fiber orientation calculated using the statistical scheme for a simple shear flow were in good agreement with known analytical results (Akbar and Altan (1992)), indicating the accuracy of the statistical scheme. In this section the number of fibers was $N = 3,600$ with aspect ratio $\lambda = 5$ and 50, and the initial orientation was isotropic.

### 4.3.2 Performance of a quadratic closure approximation

The evolution of the orientation distribution function and the orientation ellipses for $\lambda = 5$ along $\Psi = -0.9$ streamline in a developing flow are shown in Figure 19 in order to easily understand the relation between the orientation distribution function and orientation ellipse. In the computations of a developing flow of Newtonian fluids the streamfunction $\Psi$ was set equal to 0 on the centerline and -1.0 at the lower wall. In Figure 19(b) isotropic orientation (drawn by a circle) at the inlet changes into highly aligned state as the fibers flow downstream, and almost complete alignment (unit straight line) parallel to the channel wall can be seen at $t = 4.0$. Then the preferred angle of fibers becomes negative, and the orientation ellipse begins to flip over. The flip-over phenomenon of the orientation ellipse can be often observed for small aspect-ratio fibers.

Next the performance of a quadratic closure approximation is examined by comparison of **a** calculated using the statistical scheme with that obtained by the evolution equation of **a** with a quadratic closure equation. In Figure 20 solid lines with symbol represent the results computed from the statistical scheme and broken lines from the evolution equation of **a**. Fluid velocity has relatively large width-direction (y-direction) component near the channel wall within the inlet region, thus fibers firstly begin to align in the y-direction and the value of $a_{11}$ decreases from 0.5 in Figure 20(b). Then alignment in the direction parallel to the wall (x-direction) becomes gradually dominant, and the value of $a_{11}$ starts to increase. This trend becomes more pronounced as the channel wall is approached. On the other hand, the alignment in the x-direction becomes better within a region near the
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Evolution of the orientation distribution function of fibers and the corresponding orientation ellipse $\lambda = 5$ for along $\psi = -0.9$ streamline in a developing Newtonian flow; (a) orientation distribution function; (b) orientation ellipse</figcaption></figure>

centerline owing to the acceleration of fluid velocity in the *x*-direction in Figure 20(a). The above trend in the variation of **a** is in good agreement with the statistical scheme and the evolution equation of **a**. However, the quantitative agreement of the variation of **a** becomes limited. Predictions of fiber orientation in a developing flow through a parallel plate channel may clearly indicate that the much simpler quadratic closure approximation performs well when the fibers are highly aligned, however it is neither suitable for almost isotropic orientation state nor for a rapid change in fiber orientation which is caused by a significant change in fluid velocity observed within the inlet region in a developing flow or flip-over phenomenon.

The variations of the fiber orientation in a developing flow regime can be clearly seen in Figure 21: the orientation ellipse begins to align in the *x*-direction after counterclockwise rotation at the immediate downstream region of the inlet near the channel wall, while the orientation ellipse aligns well along the centerline. Furthermore, Figure 21 can exhibit the performance of a quadratic closure approximation.
<figure><img src="image_21.png" /><figcaption>**Figure 20.** Comparison between the components of the second-order orientation tensor $a_{11}$, $a_{12}$ computed from the evolution equation of **a** with a quadratic closure approximation (broken lines) and those computed by a statistical scheme (solid lines with symbols) for $\lambda = 50$ along (a) $\psi = -0.1$; (b) $\psi = -0.9$ streamlines in a developing Newtonian flow</figcaption></figure>

### 4.3.3 Conclusions

In this section the planar orientation of fibers in a developing flow of a Newtonian fluid through a parallel plate channel was analyzed in order to examine the performance of a quadratic closure approximation, and comparison of the variations of **a** calculated from the statistical scheme with the solutions of the evolution equation of **a** with a quadratic closure approximation was also presented.
<figure><img src="image_22.png" /><figcaption>**Figure 21.** Evolution of orientation ellipse in a developing Newtonian flow for $\lambda = 50$. (a) orientation ellipses computed from the evolution equation of **a** with a quadratic closure approximation; (b) orientation ellipses computed by a statistical scheme</figcaption></figure>

A quadratic closure approximation performs well when the fibers are highly aligned, however it is neither suitable for disturbed orientation state nor for a rapid change in fiber orientation which is caused by a significant change in fluid velocity or flip-over phenomenon. These results were confirmed in a simple shear flow using the solutions of the orientation distribution function by Advani and Tucker (1987), (1990). Furthermore, it was found that the statistical scheme is very useful and feasible method to analyze the orientation distribution of fibers in fiber suspension flow.

## 4.4 Analysis of Fiber Orientation in a Recirculating Flow

### 4.4.1 Computations of fiber orientation

The characterization of the fiber orientation in a recirculating flow is very important from a practical point of view. Therefore, we investigated in this work the planar fiber orientation in a recirculating flow within a salient corner of a 1:4 backward-facing step channel using both numerical simulations and fiber orientation observations. The fundamental and important results were obtained, which are necessary to rigorously compute fiber suspension flows through a complex channel.

We used the statistical scheme in the computations of fiber orientation: orientations of $N = 1800$ fibers with aspect ratio $\lambda = 5$ and $10000$ ($k \approx 1$), could be evaluated from an isotropic initial orientation along streamlines as described in subsection 4.3.

Computations of fiber orientation in a recirculating flow were continued until a fluid element containing many fibers made about 10 rounds along a streamline in order to examine the existence of the steady orientation state. In the present work the Reynolds number was defined in the downstream channel as $Re = \rho UH/\eta$.
### 4.4.2 Observations of fiber orientation

Orientation visualization experiments of large aspect-ratio fibers were carried out in Newtonian flows. Figure 22 is a schematic diagram of a backward-facing step channel used in the experiments.

Vinylon fibers of $d = 10.6~\mu\text{m}$ in diameter, $l = 3~\text{mm}$ in length ($\lambda = 283$) were suspended in 63 wt% corn syrup/water solution with volume fraction of $\phi = 2 \times 10^{-5}$. The viscosity of the Newtonian solvent is $\eta = 0.0463~\text{Pa} \cdot \text{s}$ at $16.8^{\circ}\text{C}$. In this case the flow fields were almost the same as those for the solvent because the fiber parameter $\phi\mu/\eta = 0.28$ was very small. The observation plane was the middle of the channel height and a planar sheet of light illuminated the fibers in this plane.

<figure><img src="image_23.png" /><figcaption>**Figure 22.** Schematic diagram of the fiber orientation visualization setup</figcaption></figure>

### 4.4.3 Fiber orientation in a recirculating flow

Figures 23 and 24 show fiber orientation state along three streamlines in a recirculating flow at $Re = 15$ for $\lambda = 10000$ and 5, respectively, using the orientation ellipse. Figures 23(a) and 24(a) shows the orientation state at the first circulation and Figures 23(b) and 24(b) the steady orientation state. It can be seen clearly that fiber orientation changes from random orientation (circle) at the initial position to complete orientation (straight line) as fibers move along the streamline: steady and complete orientation state can be achieved for both aspect ratios. Fibers are subjected to periodic flow field during circulation along the streamlines and they eventually reach complete orientation state, which means that all fibers in a fluid element align completely in one direction.

The effect of the aspect ratio is as follows: for large aspect-ratio fibers, all fibers align completely along the streamlines (co-linear alignment) in Figure 23(b). This result can be also observed in the experiments as shown in Figure 25. In this figure both fibers, which were obtained from two photographs at steady orientation state, and computed streamlines in a Newtonian flow are drawn. On the other hand, the preferred angle of fibers lies obliquely to the streamlines for small aspect-ratio fibers in Figure 24(b). This fiber orientation arises from flip-over phenomenon.
<figure><img src="image_24.png" /><figcaption>**Figure 23.** Evolution of orientation ellipse along $\psi = 1.01, 1.05, 1.08$ streamlines in a recirculating flow within a salient corner of a 1:4 backward-facing step channel at $Re = 15$, for $\lambda = 10000$. (a) at the first circulation; (b) at the steady orientation state. The value of the streamfunction on the vortex boundary is $\psi = 1.0$</figcaption></figure>

The degree of orientation in the preferred direction can be indicated by the maximum eigenvalue of **a**. Figure 26 shows the evolution of the maximum eigenvalue for small aspect-ratio fiber: the value of 0.5 expresses random orientation and the value of unity means the complete orientation in one direction. The values along three streamlines approach unity after starting from 0.5. Furthermore, it takes longer to reach the value of unity in the region near the vortex center. We can also confirm in the experiments that alignment of fibers with the streamlines becomes a little worse in the inner region near the vortex center (Figure 25).

In Figure 26 an oscillating orientation order parameter (maximum eigenvalue) can be clearly observed along $\psi^* = 1.08$ streamline for $\lambda = 5$. This is due to flip-over of individual fibers in a small fluid element because this streamline has a small curvature radius, then fibers tend to flip over. This phenomenon could be also observed for large aspect-ratio fibers. On the other hand, a dip in orientation order parameter which is seen along $\psi^* = 1.01$ streamline approximately at $t^* = 7$ in Figure 26 is caused by the flip-over which occurs along the vertical channel wall near the re-entrant corner shown in Figure 24(a).
<figure><img src="image_25.png" /><figcaption>**Figure 24.** Evolution of orientation ellipse along $\psi = 1.01, 1.05, 1.08$ streamlines in a recirculating flow within a salient corner of a 1:4 backward-facing step channel at $Re = 15$, for $\lambda = 5$. (a) at the first circulation; (b) at the steady orientation state. The value of the streamfunction on the vortex boundary is $\psi = 1.0$</figcaption></figure>

### 4.4.4 Conclusions

The following conclusion could be obtained from the computations and observations of fiber orientation: in a recirculating flow, for large aspect-ratio fibers both experiments and computations indicated that all fibers aligned completely along the streamlines. While, for small aspect-ratio fibers the computations predicted that complete alignment could be achieved, but the preferred angle lay obliquely to the streamlines.

Consequently, when we compute rigorously complex flows of fiber suspensions by coupling flow field with fiber orientation, co-linear alignment of fibers is valid for large aspect-ratio fibers in a recirculating flow. While the preferred angles do not coincide with the streamlines in a recirculating flow for small aspect-ratio fibers, so these angles have to be determined from the steady orientation of a single fiber. In this case we do not need to use the statistical scheme because of complete alignment of fibers, then the computational efficiency is greatly enhanced.
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Observed fiber orientation state and computed streamlines in a recirculating flow within a salient corner of a 1:4 backward-facing step channel at $Re = 14.9$</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 26.** Evolution of orientation order parameter in a recirculating flow within a salient corner of a 1:4 backward-facing step channel at $Re = 15$, for $\lambda = 5$</figcaption></figure>
## 4.5 Numerical Solution of Fiber Suspension Flow Through a Rectangular Channel

The coupled flow kinematics and fiber orientation distribution were computed to study the development of fiber suspension flows through a parallel plate channel. Fibers can align in three-dimensional way even in two-dimensional flows. In this work, therefore, we computed the two-dimensional flows of suspensions with high aspect-ratio fibers in a Newtonian fluid for planar and three-dimensional fiber orientation distributions. In this subsection comparison of flow field, fiber orientation and stress distribution between for the planar orientation of fibers and for the three-dimensional case will be described.

### 4.5.1 Solution procedure

We will briefly explain the solution procedure. The continuity and momentum equations for the flow are Eq. (4.1) and Eq. (4.2), respectively, while the extra stress for the fiber suspension is determined using Eq.(4.3). All equations were written in terms of non-dimensional variables using the channel width $H$ and the mean fluid velocity $U$ as the reference parameters. Subsequently, the Reynolds number was defined as $Re = \rho U H / \eta$. The computational domain was $x^* = 0$ to 10 in the flow direction and $y^* = 0$ to 0.5 in the width direction (a lower half part of the channel).

The development of the flow field and planar/three-dimensional fiber orientation field was studied for the initial parabolic velocity profile and isotropic orientation at the channel inlet. In the present study, the statistical scheme was employed. In the computations, $k$ was set equal to unity.

The components of the fourth-order orientation tensor $a_4$ can be calculated by Eq. (4.13) for the planar orientation and by Eq. (4.14) for the three-dimensional orientation of fibers

$$ 
\begin{aligned}
a_{1111} &= \sum_{i=1}^{N} \cos^4 \beta_i / N, & a_{1112} &= \sum_{i=1}^{N} \cos^3 \beta_i \sin \beta_i / N \\
a_{1122} &= \sum_{i=1}^{N} \cos^2 \beta_i \sin^2 \beta_i / N, & a_{1122} &= \sum_{i=1}^{N} \cos \beta_i \sin^3 \beta_i / N \\
a_{2222} &= \sum_{i=1}^{N} \sin^4 \beta_i / N
\end{aligned} 
\quad (4.13) $$

and

$$ 
\begin{aligned}
a_{1111} &= \sum_{i=1}^{N} \sin^4 \theta_i \cos^4 \beta_i / N, & a_{1112} &= \sum_{i=1}^{N} \sin^4 \theta_i \cos^3 \beta_i \sin \beta_i / N, \\
a_{1122} &= \sum_{i=1}^{N} \sin^4 \theta_i \cos^2 \beta_i \sin^2 \beta_i / N, & a_{1122} &= \sum_{i=1}^{N} \sin^4 \theta_i \cos \beta_i \sin^3 \beta_i / N \\
a_{2222} &= \sum_{i=1}^{N} \sin^4 \theta_i \sin^4 \beta_i / N
\end{aligned} 
\quad (4.14) $$

where $\beta$ is an angle between the $x$-axis and the projection of the fiber in the flow plane ($x-y$ plane), and $\theta$ is an angle between the $z$-axis and the fiber.

As already described in section 4.3, the number of fibers of 1800 or more was appropriate for planar orientation. However, in the present computations the number of fibers was
<figure><img src="image_28.png" /><figcaption>**Figure 27.** Effect of number of fibers on the component of the fourth-order orientation tensor $a_{1111}$ for the three-dimensional orientation of fibers in a developing flow of fiber suspension through a parallel plate channel at $Re = 10$, for $\phi\mu/\eta = 10$</figcaption></figure>

$N = 180$ because of the limitation of computational time. On the other hand Figure 27 shows the effect of number of fibers on $a_{1111}$ for the three-dimensional orientation of fibers: the number of segmentations between 0 and 90° in the $\theta$ direction were varied from $N_s = 9$ to $N_s = 90$. In this case, the total number of fibers was varied from 202 to 20,582. The values of $a_{1111}$ for $N_s = 30$ to 90 are in good agreement. Then, $N_s$ was set equal to 30 ($N = 2277$) in the computations.

### 4.5.2 Velocity profile and fiber orientation

Comparison of centerline velocity profile between the planar orientation and the three-dimensional orientation case is shown in Figure 28. The centerline velocity profile for fiber suspensions can change most significantly from the Newtonian case at the location of approximately $x^* = 0.9 \sim 0.975$. However, the effect of the addition of short fibers to a Newtonian solvent on the centerline velocity profile gradually disappears as the downstream region is approached. Figure 28 shows clearly that the present computational domain, length of 10, is not long enough for the flow to develop fully. However, we can easily suppose that the flow kinematics for fiber suspensions finally become identical to the Newtonian flow in a channel long enough to reach a fully-developed flow. This phenomenon is due to complete alignment of fibers parallel to the channel wall in the fully-developed flow regime.

The degree of isotropy of fiber orientation at the inlet is stronger for the three-dimensional orientation case than for planar one, therefore, Figure 28 shows clearly that the deviations of the centerline velocity profile for fiber suspension from Newtonian counterpart become
<figure><img src="image_29.png" /><figcaption>**Figure 28.** Centerline velocity profiles in a developing flow of fiber suspension through a parallel plate channel at $Re = 10$; for the planar orientation (solid lines); for the three-dimensional orientation (thick dotted lines)</figcaption></figure>

less significant for the three-dimensional orientation of fibers. However, the entrance region necessary for the flow to fully develop becomes longer for the three-dimensional orientation of fibers.

This result is confirmed by comparing of the evolution of $a_{1111}$ between the three-dimensional orientation and the planar orientation of fibers (Figure 29). The value of $a_{1111}$ becomes smaller for the three-dimensional orientation case.

Figure 29 shows the variations of $a_{1111}$ in the width direction. At the inlet the orientation of fibers is isotropic except parallel alignment at the channel wall, thus $a_{1111} = 3/8$ for the planar orientation and $a_{1111} = 1/5$ for the three-dimensional orientation. Fibers align gradually in the x direction as they flow down to the exit of the channel. Therefore, $a_{1111}$ approaches unity within almost all region near the exit. However, a long distance is necessary for fibers to change their orientations from isotropic distribution to complete alignment in the region near the centerline owing to low velocity gradient.

The effects of the fiber parameter on the profiles of the x component of velocity $u^*$ in the width direction can be seen in Figure 30. As the addition of short fibers to Newtonian liquids, the velocity profile changes from the parabolic profile to a more plug-like one because fibers tend to align in the x direction and the flow drag is reduced as the wall is approached. For the planar orientation, the velocity profile becomes more plug-like than that for the three-dimensional orientation of fibers, because the degree of isotropy for the planar orientation is weaker than that for the three-dimensional orientation case.

### 4.5.3 Stress distribution in fiber suspension flow

Finally, the shear stress $T_{xy}^*$ and the normal stress difference $T_{xx}^* - T_{yy}^*$ are illustrated within $0 \le x^* \le 2$ in Figures 31 and 32, respectively. The non-dimensional stress is defined as
<figure><img src="image_30.png" /><figcaption>**Figure 30.** Velocity profiles in the width direction in a developing flow of fiber suspension through a parallel plate channel at $Re = 10$, for $\phi\mu/\eta = 10$; for the planar orientation (solid lines with symbols and thick broken lines); for the three-dimensional orientation (color lines) at $x^* = 0, 0.5, 1.25, 2.5, 10$</figcaption></figure>

The anisotropic characteristics of fiber orientations and stress field for fiber suspensions were remarkably observed in the region near the inlet. As a result, the flow kinematics and stress field for fiber suspensions changed more significantly from the Newtonian counterpart in the region near the inlet. These results arise from the fact that a parabolic velocity profile was given at the inlet of the channel. However, the effect of the addition of short fibers gradually disappeared as the downstream region was approached.

The degree of isotropy for the three-dimensional orientation is stronger than that for the planar orientation case near the inlet, then, the velocity profile in the width direction became more plug-like and also the variations of the stresses were more significant for the planar orientation of fibers.

## 4.6 Concluding Remarks

Fiber suspension flow and fiber orientation through a complex channel on the basis of both numerical simulation and visualizations of flow and fiber orientation were described in this section. We focused on the case of rigid fiber suspensions in Newtonian liquids, and examined four particular topics: (1) Entry flow of dilute fiber suspensions through a circular contraction; (2) Discussion on a quadratic closure approximation using a statistical scheme; (3) Analysis of fiber orientation in a recirculating flow; (4) Numerical solution of fiber suspension flow through a rectangular channel.

The most obvious result in an entry flow of fiber suspensions is that the salient corner vortex significantly grows and also the normal stress difference largely increases for fiber suspensions as compared to Newtonian fluid even at a very low volume fraction of fibers.

A statistical scheme was confirmed to be a very useful and feasible method to analyze
<figure><img src="image_31.png" /><figcaption>**Figure 31.** Distribution of the shear stress $T_{xy}^*$ in a developing flow of fiber suspension through a parallel plate channel at $Re = 10$, for $\phi\mu/\eta = 10$; (a) for the planar orientation, the maximum value and the location are 14.4 (0.025, 0.025); (b) for the three-dimensional orientation, 11.1 (0.025, 0.025)</figcaption></figure>

accurately the orientation distribution in fiber suspension flow. In the statistical scheme the orientations of large numbers of fibers are evaluated from computation of the Jeffery equation along the streamlines instead of direct solutions of the orientation distribution function of fibers or the evolution equation of the orientation tensor with a closure equation. We, therefore, used the statistical scheme to examine the performance of a quadratic closure approximation, investigated the fiber orientations in a recirculating flow and computed fiber suspension flows through a rectangular channel by coupling flow kinematics and fiber orientation distribution.

1. A quadratic closure approximation performs well when the fibers are highly aligned, however it is neither suitable for disturbed orientation state nor for a rapid change in fiber orientation.

2. In a recirculating flow, all fibers align completely along the streamlines (co-linear alignment) for large aspect-ratio fibers, while complete alignment can be achieved but the preferred angle lies obliquely to the streamlines for small aspect-ratio case.

3. The flow kinematics and stress field for fiber suspensions in a developing flow through a parallel plate channel change more significantly from the Newtonian counterpart in the region near the inlet. These results arise from a parabolic velocity profile specified at the inlet of the channel.

On the other hand, the future works on fiber suspension flows are as follows: the inter-
<figure><img src="image_32.png" /><figcaption>**Figure 32.** Distribution of the normal stress difference $T_{xx}^* - T_{yy}^*$ in a developing flow of fiber suspension through a parallel plate channel at $Re = 10$, for $\phi\mu/\eta = 10$; (a) for the planar orientation, the maximum and minimum values and the location are 13.9 (0.075, 0.025), -3.44 (0.025, 0.5); (b) for the three-dimensional orientation, 10.2 (0.075, 0.025), -1.82 (0.025, 0.5)</figcaption></figure>

actions of both fiber-fiber and fiber-wall are key factors for concentrated fiber suspensions. We, therefore, have to evaluate the interactions for future computations of fiber suspension flows. Furthermore, from a practical point of view the importance of fiber orientations in viscoelastic flows has been recognized. Ait-Kadi and Grmela (1994) and Azaiez (1996) proposed the constitutive equation for fiber suspension in polymeric liquids and computed steady and transient shear flow properties. Contribution of the interaction between fibers and polymer molecules to the stress of fiber suspension was taken into account in their constitutive equation, however effect of polymer molecules on the evolution of fiber orientation remains an open question.

# 5 EFFECTS OF FIBER ADDITIVES ON FLOW INSTABILITIES AND THE DEVELOPMENT OF TURBULENCE

It is well known that spectacular reduction in turbulent friction can be obtained through the addition of small amounts of polymers, polyelectrolytes, surfactants or rigid particles. This phenomenon, commonly termed as drag reduction, is a resurging subject of research which is gaining a lot of interest. The early reviews by Hoyt (1972) and Berman (1978) have documented a large number of experimental observations that bear evidence to the great potential of such additives in reducing turbulent effects.

Most of the existing studies are devoted to polymer additives that have been the subject of extensive theoretical and experimental investigations (See Lumley 1969 and the recent
review by Tiederman 1990). However, in spite of the large amount of studies dealing with polymer flows, the exact mechanism of drag reduction is still unclear. In particular, the role of stress anisotropy due to polymer extension versus elasticity, in the mechanism of drag reduction is still an ongoing subject of controversy. De Gennes (1990) and Joseph (1990), among other authors, suggest that elasticity is the main factor behind drag reduction by polymer additives. Other authors maintain that the main mechanism for drag reduction in polymers is associated with the anisotropy arising from the particular orientation of the polymer chains once they are fully stretched (Landahl 1973, Hinch 1977 and Draad and Hulsen 1995). Experiments by Sasaki (1991a,b) involving the use of polymers in combination with different types of solvents suggest that the anisotropy associated with a rod-like behaviour of the additives is essential for the drag reduction.

Studies of the effects of fiber additives on the stability and the development of turbulent flows are rather limited. Unlike their polymer counterpart, fiber additives have received very limited attention, and most of the early studies have been carried out to address specific problems encountered in practical industrial processes, and therefore were published in specialized journals associated with certain industries such as the pulp paper or textiles. However, these early studies showed clearly that fiber additives can induce drag reduction effects that are as spectacular as those obtained by polymer additives. Furthermore, it became apparent from these studies that fiber additives present a number of practical advantages over polymer additives. Indeed, fibers have the advantage of higher resistance to chemical and mechanical degradation that makes them more adapted to practical applications. Furthermore, fiber additives are more easily recoverable from the suspending fluid and their size and shape are more easily controllable. These practical advantages of fiber additives, added to their strong potential for reducing flow resistance, make them an excellent drag reducing agent.

More recently, there seems to be a resurgent interest in studying flows of fiber suspensions and their role in drag reduction. This interest can be attributed in part to recent suggestions that polymer and fiber additives affect the flow turbulence in the same way. Indeed, few similarities reported in some experimental studies involving polymer and fiber additives, and the fact that polymer chains undergoing strong strains may act hydrodynamically like solid particles, has led many authors to speculate that these two additives modify flow instabilities in a similar manner (Landahl 1973, Virk and Wagger 1990, and Draad and Hulsen 1995). Unfortunately there is still no clear proof that the effects of these two additives on flow turbulence are the same, and it is suspected that they may involve fundamentally different mechanisms. Indeed there are many reasons why the analogies between polymer and fiber additives should not be taken too far. First, rigid particles lack the flexibility that polymer chains usually have and which may have important hydrodynamic effects. Second, it is not obvious that a stretched polymer chain can retain its conformation without degradation, and therefore behave exactly like a long rigid particle. Furthermore, it is known that in general, for a given additive concentration, polymers tend to induce stronger drag reduction effects than their solid particle counterparts. These differences are confirmed by drag reduction measurements reported by Lee, Vaseleski and Metzner (1974) for pipe flows, that show that the combined use of polymers and fibers leads to much larger drag reduction effects than what is usually obtained with either additive alone. Furthermore, laser-Doppler anemometry measurements by McComb and Chan (1979) indicate that effects of fibers on the turbulent structures in a pipe flow were fundamentally different from those of polymer additives.

This section of the present review is devoted to examining the role of fiber additives in changing the mechanisms of instability and development of the full turbulent regime in a wide variety of flows encountered in nature and industry. Because of the brevity of this section, the literature cited is subjectively limited though it covers most relevant studies.
We will first review studies dealing with pipe flows. Then we will focus on the effect of fiber additives on the centrifugal instability in Taylor-Couette flows. Finally, we will examine the special case of free shear flows.

## 5.1 Pipe Flows

It has been known that the presence of fibers in pipe and channel flows can result in a large reduction of the frictional drag. This effect is usually measured by the drag reduction defined as:

$$Dr = \frac{f_w - f_s}{f_w} \qquad (5.1)$$

where $f_w$ and $f_s$ are the friction factors of the suspending fluid in the absence of fibers and the suspension, respectively. The friction factor is defined as

$$f = \frac{\tau}{\frac{\rho v^2}{2}} \qquad (5.2)$$

where $\tau$ is the shear stress, $v$ the bulk flow velocity and $\rho$ the density of the suspending fluid.

As mentioned earlier, most of the existing literature dealing with drag reduction by additives has mainly focused on the use of polymer additives particularly in wall-bounded flows, due to their importance in many technological processes (see for example Pinho and Whitelaw 1990, and Tiederman 1990). On the other hand, the use of fiber additives as potential drag reducing agents remains very limited. Nevertheless, ancient and recent experimental studies show conclusively that fibers can be a very effective drag reducing agent in a wide variety of flows.

Among the earliest studies dealing with the flow of fiber suspension, we will mention the works of Daily and Bugliarello (1961) and Arranaga (1970) and Ellis (1970). These experimental studies allowed to examine the effects of fiber additives on the flow dynamics in pipe flows. The characterization of the flow resistance revealed that the friction factor is a monotonically decreasing function of the Reynolds number determined using the viscosity of the suspending fluid (water). Furthermore, the authors identified three flow regimes in which the fibers have a distinctive effect on the flow. In the first regime that develops at low flow rates, the flow consists of a central core where most of the fibers are located, surrounded by an annulus consisting of an almost fiber-free suspending fluid. The examination of the velocity profiles shows that most of the velocity gradient occurs in the peripheral annulus. In this regime, the presence of fibers in the flow results in stronger flow resistance, and the friction factor increases with the fiber concentration and aspect ratio. In the transition regime, fibers start to be separated from the plug flow in the central core and the velocity gradients near the wall become less steep. In this regime the friction factor can be lower than that for pure water. Finally, in the turbulent regime the plug is broken and fibers are spread throughout the pipe. The friction coefficient is less than that for the pure carrier fluid, and varies little with the Reynolds number.

Mean velocity profile measurements were conducted by Mih and Parker (1967), who noticed that a plug core flow forms in the center of the pipe due to the presence of fibers. It was also observed that the diameter of the plug core decreases with decreasing fiber concentration and increasing velocity.

Pirih and Swanson (1972) used hot film anemometry to study the effects of small rigid particles on the turbulent structures in pipe flows. It was reported that unlike pure water flows which exhibit an abrupt increase in the friction factor as the flow becomes turbulent, suspension flows resulted in a smooth transitional behavior with a continuously decreasing
friction factor. The authors studied also the turbulence intensity and compared the longi-
tudinal energy spectra for fiber suspensions and pure water in pipe flows. It was found that
for the same flow velocity, suspension flows develop larger eddies and more stable turbulent
structures than pure water flows.

Vaseleski & Metzner (1974) conducted a series of experimental studies to determine the effects of fiber additives on turbulent pipe flows. They found that turbulent friction effects were reduced as the fiber aspect-ratio or volume fraction is increased. A subsequent study by Kale and Metzner (1976) showed that drag reduction due to fiber additives is a wall region phenomenon. The effects of the fibers are confined to the region of the flow close to the pipe wall where eddies are mainly of dissipative type and form a viscous zone of flow turbulence. This result was confirmed by Sharma *et al.* (1978) who studied drag reduction effects resulting from the injection of fibers in the flow. The authors used two modes of injection where the fibers where either injected close to the pipe wall or near the centerline. It was found that the boundary injection leads to stronger drag reduction effects in turbulent pipe flows than the centerline injection. This study also suggests that adding fibers through injection into the flow leads to stronger drag reduction effects than the traditional approach using premixed suspensions, and that using the injection technique, strong drag reduction effects can be obtained even with a trace quantity of fibers.

A study by McComb and Chan (1979) who investigated a pipe flow of a 300 wppm suspension of high aspect ratio asbestos fibers, allowed to measure the streamwise velocity profiles over cross sections of the pipe using LDA. It was reported that due to fibers degradation, the drag reduction effects decreased and reached levels comparable to those of polymer additives. It was concluded that due to degradation, there is a transition from *fiber-like* effect where the drag reduction effects are strong and the streamwise intensity are below that of the solute alone to a *polymer-like* behaviour with increased streamwise intensity and weaker drag reduction. These results were further expanded in a study that presents measurements of the streamwise and circumferential velocities as well as the energy spectra of the flow (McComb and Chan 1985).

5.2 Taylor Flow

A Taylor-Couette flow develops between two coaxial circular cylinders of radii $R_1$ and $R_2$ ($R_2 > R_1$), rotating about the common axis with angular velocities $\omega_1$ and $\omega_2$. This type of flow is of special interest from both fundamental and practical points of view.

In the case where the outer cylinder is fixed, a circular Couette flow develops at low Reynolds number. As $Re$ is increased, a transition is observed at a critical Rec wherein toroidal counter-rotating cells separated by transversal planes develop. This flow is known as the Taylor vortex flow. At an even higher critical value of the Reynolds number, the horizontal cells become wavy, this is the wavy Taylor vortex flow. Other flow transitions can be observed at higher Reynolds number, such as the modulated wavy Taylor vortices, but the flow quickly becomes turbulent. It is worth noting that the above transitions can be also observed under a wide range of conditions when both cylinders are rotating whether in the same sense or in opposite senses.

The first attempts to analyze the stability of a Taylor-Couette flow in the presence
of particle additives used a theory of incompressible anisotropic fluids proposed by Erick-
sen (1960) to describe the flow. The objective of these studies was to determine whether
the presence of fiber additives reduces or enhances the instability of the flow. The main
weakness of these studies lays in the constitutive model used to determine the fiber orien-
tation and its contribution to the total stress. Indeed, due to the mathematical complexity
in describing the fiber orientation, most studies have used the fiber-aligned assumption
to determine the local fiber orientation. In this approach, it is assumed that the fibers
are everywhere aligned with the streamlines such that the fiber contribution to the bulk
stress is directly expressed in terms of the velocity field instead of solving explicitly for the orientation distribution. However, because of the inevitable presence of hydrodynamic interactions, it is a truism to say that the orientation of the fibers will not be necessarily along the streamlines everywhere in the flow.

The earliest studies of the stability of the circular Couette flow are attributed to Leslie (1964) and Ericksen (1966), who examined the flow using the small gap limit where $R_1/R_2 \approx 1$ and focused on non-Brownian particles.

Pilipenko et al. (1981) used the Ericksen model to investigate the flow stability of a suspension of weakly Brownian slender fibers in the gap between two coaxial cylinders. These authors used linear stability analysis and experimental measurements to determine the drag reduction effects due to fiber additives. It was found that in the small gap limit, the addition of fibers to the flow has a destabilizing effect and leads to a decrease in the critical Taylor number at which the secondary flow develops.

Nsom (1994-a) conducted a linear stability analysis of the Taylor-Couette flow occurring between two concentric cylinders in the presence of stiff macroscopic fibers. Unlike previous studies that used the small gap approximation, this analysis considered the situation of wide-gap geometry to account for the macroscopic nature of the fibers, and examined the dilute and semi-concentrated regimes. For the latter regime, the author modified the Ericksen model using the formulation proposed by Shaqfeh et al. (1990) to account for fiber-fiber interactions. It was found that the presence of fibers does not affect the critical wave-number, but does lead to larger critical Taylor numbers when the fiber aspect-ratio or volume fraction are increased. In a subsequent study, Nsom (1994-b) computed the drag reduction for the Taylor-vortex flow for dilute and semi-dilute fiber suspensions. It was found that for dilute suspensions, the drag reduction increases with the aspect ratio for a given cylinder radii ratio. In the case of semi-dilute suspensions, it was reported that the drag reduction increases with the volume fraction for a given cylinder radii ratio and fiber density.

A similar study was conducted for the Dean flow occurring in a curved channel under a pressure gradient (Nsom, 1996). The results of the analysis were similar to those obtained for the Taylor-Couette flow, suggesting that an increase in the fiber aspect-ratio or volume fraction tends to reduce the flow instability.

A more recent study by Gupta et al. (2001) investigated the linear stability of rigid non-Brownian fibers in a Taylor-Couette flow. These authors determined explicitly the fiber orientation by using a constitutive model based on orientation tensors in conjunction with the hybrid closure. It was found that the presence of fibers in the flow reduces the flow instability, and that the critical Reynolds number increases with increasing inter-fiber interaction coefficient, volume fraction or aspect ratio. This suppression of the flow instability was attributed to the development of negative first normal stress difference that result in a state of compression of fluid elements in the radial direction. Furthermore, an analysis of the budget of the disturbance kinetic energy revealed that the orientation perturbation plays a dissipative role in the flow.

## 5.3 Free Shear Flows

Free shear flows such as mixing layers, jets and wakes differ from the previously mentioned wall bounded flows in that they have an inflectional mean velocity profile, and hence are subject to inviscid instabilities. These flows are encountered in a wide variety of natural and technological systems and it is important to understand the mechanisms governing the process of transition to turbulence in order to predict, and if possible control the evolution of such flows.

All existing experimental investigations of free shear flows with additives, whether viscoelastic substances or rigid particles, bear evidence to the great potential of such additives
to reduce these flows instability. However, in spite of the overwhelming experimental and theoretical evidence for drag reduction by such additives, the physical mechanisms responsible for this phenomenon are not completely understood and remain a subject of debate.

The stability of this class of flows in the presence of fiber additives has been the subject of a very limited number of studies. The experimental visualizations of a jet flow by Filipsson et al. (1977) represent one of the few experiments devoted to the study of the stability of free shear flows in the presence of particle additives. The authors reported that the addition of a small amount of fibers (asbestos) led to drastic changes in the flow that translated in an enhancement of large-scale turbulent structures and a modulation of the turbulence by the suppression of small-scale structures. It was also observed that similar trends were obtained when a polymer (Polyox WSR-301) was added to the flow, which led the authors to speculate that the polymer and fiber additives modify the Kelvin-Helmholtz instability in a similar manner.

Yang et al. (1990) examined the spatial and temporal stability of a gaseous mixing layer in the presence of rigid spherical particles. It was found that the presence of particles decreases the amplification rate of the perturbation therefore reducing the instability of the flow. It was reported that the reduction of the flow instability is stronger for higher particle loading, and that the maximum growth rate decreased almost linearly with the particle loading.

Recently, Azaiez (2000a) conducted temporal linear stability analyses of a liquid mixing layer flow in the presence of elongated rigid particles. The author used a formulation based on moments of the probability distribution function to determine the particle orientation, and extended the classical linear stability theory to derive modified Orr-Sommerfeld equation. To determine the expression of the fourth-order orientation tensor, hybrid and quadratic closure approximations for 2D orientation were used:

$$ \textit{Quadratic:} \quad a_{ijkl} = a_{ij}a_{kl} $$

$$ \textit{Linear:} \quad a_{ijkl} = \frac{\det(\mathbf{a})}{6}(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \frac{1}{3}(a_{ij}a_{kl} + a_{ik}a_{jl} + a_{il}a_{jk}) $$

$$ \textit{Hibrid:} \quad a_{ijkl} = f a_{ij} a_{kl} + (1 - f) b_{ijkl} \qquad (5.3) $$

$$ b_{ijkl} = -\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}) $$

$$ f = 1 - 4 \det(\mathbf{a}) $$

The linear stability analysis used a tanh base state velocity profile, and the components of the orientation tensors were determined accordingly (Azaiez, 2000-a). The temporal stability of the flow was investigated using the classical normal mode analysis (Drazin and Reid, 1987): The base flow is perturbed by infinitesimal disturbances, whose evolution is examined by linearizing the governing equation about the base state. The perturbations are decomposed into Fourier modes as follows:

$$ \omega(x, y, t) = \omega_0(y) + \omega'(x, y, t) = \omega_0(y) + \omega^*(y)e^{i(ax - \sigma t)} \qquad (5.4) $$

$$ \mathbf{a}(x, y, t) = \mathbf{a}_0 + \mathbf{a}'(x, y, t) = \mathbf{a}_0 + \mathbf{a}^*(y)e^{i(ax - \sigma t)} $$
If the disturbance is found to grow with time $Im(\sigma) > 0$, the flow is considered stable while the flow is regarded as stable if the perturbation is damped $Im(\sigma) < 0$. A zero growth rate, $Im(\sigma) = 0$, corresponds to a neutrally stable flow. This approach allows to transform the original set of partial differential equations into an ordinary differential eigenvalue problem. Three different methods were used to solve the linearized problem with the appropriate boundary conditions. The first one consisted of an iterative method based on orthogonal shooting which has been described in details in the paper of Azaiez and Homsy (1994). The second method is based on a multiple shooting technique (Ascher et al. 1988). In this method, one divides the whole interval by a mesh into $N$ sub-intervals and shoots in each interval. The final result is obtained from the different solutions in the sub-intervals by imposing matching conditions at the shooting points. This method is very useful when the problem has both rapidly decaying and rapidly growing solutions. In the third method, we solved the linearized vorticity and fiber orientation equations using the finite difference method. The ordinary differential equations are transformed into a set of finite difference equations and lead to a generalized eigenvalue problem $A \cdot \vec{x} = \lambda B \vec{x}$, which was solved using a standard QR algorithm (Azaiez 2000-a). For all three methods, the accuracy of the results was tested by refining the mesh and varying the width of the domain in the transverse direction.

This study revealed that an increase in the fiber aspect-ratio, volume fraction or hydrodynamic interaction coefficient leads to a strong attenuation of the flow instability. Both the maximum growth rate of the disturbance and the spectrum of unstable wavenumbers are reduced. It also showed that the stabilizing effects arise from the orientation diffusion due to hydrodynamic interactions, and not from the anisotropy induced by the presence of fibers in the flow, as speculated before. An examination of the budget equation for the rate of production of entropy, revealed that the main factor behind the reduction of the flow instability is associated with the fiber shear stress disturbance. This disturbance is important when hydrodynamic interactions are included in the flow, and acts as a dissipative term as the fibers tend to deviate from the fully aligned anisotropic orientation. It was also concluded that fiber models based on the assumption that particles are fully aligned with the flow and therefore ignore any misalignment induced by hydrodynamic interaction, will fail to capture the main physics of the flow.

In a subsequent study, we investigated the dependence of the results of the stability analysis of the previous study on the nature of the closure approximation, as well as the effects of fiber off-plane orientation (Azaiez, 2000-b). In particular, the stability analysis was undertaken using a 3D hybrid closure approximation :

$$ 
\begin{aligned}
\text{Hibrid (3D): } & a_{ijkl} = fa_{ij}a_{kl} + (1-f)b_{ijkl} \\
& b_{ijkl} = -\frac{1}{35}(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \\
& \qquad +\frac{1}{7}(a_{ij}\delta_{kl} + a_{ik}\delta_{jl} + a_{il}\delta_{jk} + a_{kl}\delta_{ij} + a_{jl}\delta_{ik} + a_{jk}\delta_{il}) \\
& f = 1 - 27 \det(\mathbf{a})
\end{aligned}
\qquad (5.5) $$

It was found that accounting for off-plane orientation by using the 3D hybrid closure approximation leads to substantially different instability characteristics. In particular, the flow was considerably less unstable than its planar counterpart. Furthermore, the flow dynamics were fundamentally different in the 2D and 3D cases. In the case of a 2D orientation, the shear stress disturbances tend to enhance the flow stability while the normal
stress disturbances act towards increasing the instability. On the other hand, for a 3D orientation both shear and normal components of the stress disturbance act as a stabilizing factor in the flow. This fundamental difference between the 2D and 3D configurations was attributed to the differences in the resulting orientation distributions.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge the reproduction from Journal of Non-Newtonian Fluid Mechanics, Vol. 78, K. Chiba and K. Nakamura, “Numerical solution of fiber suspension flow through a complex channel”, pp. 167-185, Copyright (1998), and Vol. 99, K. Chiba, K. Yasuda and K. Nakamura, “Numerical solution of fiber suspension flow through a parallel plate channel by coupling flow field with fiber orientation distribution”, pp. 145-157, Copyright (2001), with permission from Elsevier Science. The author gratefully acknowledge R. Torres and I. Monton, who have codified the 3D simulation code used in the simulation of the short fiber composites injection processes.

## REFERENCES

1. T. Abdul-Karem, D.M. Binding and M. Sindelar (1993), “Contraction and expansion flows of non-Newtonian fluids”, *Compos. Manuf.*, **4**, 109.

2. S.G. Advani and C.L. Tucker III (1987), “The use of tensors to describe and predict fiber orientation in short fiber composites”, *J. Rheol.*, **31**, 751.

3. S.G. Advani and C.L. Tucker III (1990), “Closure approximations for three-dimensional structure tensors”, *J. Rheol.*, **34**, 367.

4. A. Ahmed and A.N. Alexandrou (1994), “Unsteady flow of semi-concentrated suspensions using finite deformation tensors”, *J. Non-Newtonian Fluid Mech.*, **55**, 115.

5. A. Ait-Kadi and M. Grmela (1994), “Modelling the rheological behaviour of fibre suspensions in viscoelastic media”, *J. Non-Newtonian Fluid Mech.*, **53**, 65.

6. S. Akbar and M.C. Altan (1992), “On the solution of fiber orientation in two-dimensional homogeneous flows”, *Polym. Eng. Sci.*, **32**, 810.

7. M.C. Altan, S.G. Advani, S.I. Guceri and R.G. Pipes (1989), “On the description of the orientation state for fiber suspensions in homogeneous flows”, *J. Rheol.*, **33**, 1129.

8. M.C. Altan, S.I. Guceri and R.B. Pipes (1992), “Anisotropic channel flow of fiber suspensions”, *J. Non-Newtonian Fluid Mech.*, **42**, 65.

9. M.C. Altan and B.N. Rao (1995), “Closed-form solution for the orientation field in a centered disk”, *J. Rheol.*, **39**, 581.

10. U.M. Ascher, R.M.M. Mattheij and R.D. Russel (1988), “Numerical solution of boundary value problems for ordinary differential equations”, Series in *Computational Mathematics*, Prentice Hall, Englewood Cliffs, N.J.

11. A.B. Arranaga (1970), “Friction reduction characteristics of fibrous and colloidal substances”, *Nature*, **225**, 447.

12. G. Ausias, J.F. Agassant and M. Vincent (1994), “Flow and fiber orientation calculations in reinforced thermoplastic extruded tubes”, *Intern. Polym. Proc.*, **9**, 51.

13. J. Azaiez and G.M. Homsy (1994), “Linear stability of free shear flow of viscoelastic liquids”, *J. Fluid Mech.*, **268**, 37.

14. J. Azaiez (1996), “Constitutive equations for fiber suspensions in viscoelastic media”, *J. Non-Newtonian Fluid Mech.*, **66**, 35.
15 J. Azaiez, R. Guenette and A. Ait-Kadi (1997), "Investigation of the abrupt contraction flow of fiber suspensions in polymeric fluids", *J. Non-Newtonian Fluid Mech.*, **73**, 289.

16 J. Azaiez (2000a), "Linear stability of free shear flows of fiber suspensions", *J. Fluid Mech.*, **404**, 179.

17 J. Azaiez (2000b), "Reduction of free shear flows instability: effects of polymer versus fiber additives", *J. Non-Newtonian Fluid Mech.*, **91**, 233.

18 A. Baloch and M.F. Webster (1995), "A computer simulation of complex flows of fiber suspensions", *Comput. Fluids*, **24**, 135.

19 G.K. Batchelor (1970a), "The stress system in suspension of force free particules", *J. Fluid Mech.*, **41**, 545.

20 G.K. Batchelor (1970b), "Slender bodies theory for particles of arbitrary cross section in Stokes flows", *J. Fluid. Mech.*, **44**, 419.

21 N. Berman (1978), "Drag reduction by polymers", *Ann. Rev. Fluid Mech.*, **10**, 47.

22 G. Chaidron and F. Chinesta (2002a), "On the steady solution of non-linear advection equations in steady recirculating flows", *Computer Meth. Appl. Mech. Eng.*, **191**, 1159.

23 G. Chaidron and F. Chinesta (2002b), "Further considerations in steady recirculating flows", *J. Non-Newtonian Fluid Mech.*, In press.

24 R.Y. Chang, F.C. Shiao and W.L. Yang (1994), "Simulation of director orientation of liquid crystalline polymers in 2-D flows", *J. Non-Newtonian Fluid Mech.*, **55**, 1.

25 K. Chiba, K. Nakamura and D.V. Boger (1990), "A numerical solution for the flow of dilute fiber suspensions through an axisymmetric contraction", *J. Non-Newtonian Fluid Mech.*, **35**, 1.

26 K. Chiba and K. Nakamura (1994), "Numerical solution of anisotropic channel flow. Part 2: Fiber orientation in a developing flow through a parallel plate channel", *Trans. J. Text. Mach. Soc. Japan*, **47**, T291-T299.

27 K. Chiba, T. Yamamoto and K. Nakamura (1994), "Entry flow of dilute fiber suspensions through a tubular contraction", *Nihon Reoroji Gakkaishi*, **22**, 98.

28 K. Chiba and K. Nakamura (1995), "On numerical solution of anisotropic channel flow. Part 3: Discussion on a closure approximation", *Nihon Reoroji Gakkaishi*, **23**, 159.

29 K. Chiba and K. Nakamura (1996a), "Numerical solution of fiber orientation in a recirculating flow. Part 1: Orientation of a single fiber in a recirculating flow within a slot", *Trans. J. Text. Mach. Soc. Japan*, **49**, T85-T91.

30 K. Chiba and K. Nakamura (1996b), "Numerical solution of fiber orientation in a recirculating flow. Part 2: Analysis of fiber orientation in a recirculating flow within a slot using a statistical method", *Trans. J. Text. Mach. Soc. Japan*, **49**, T107-T113.

31 K. Chiba, K. Kawata and K. Nakamura (1997), "Numerical solution of fiber orientation in a recirculating flow. Part 3: Fiber orientation in a recirculating flow within a salient corner of an abrupt expansion channel", *Trans. J. Text. Mach. Soc. Japan*, **50**, T1-T8.

32 K. Chiba and K. Nakamura (1998), "Numerical solution of fiber suspension flow through a complex channel", *J. Non-Newtonian Fluid Mech.*, **78**, 167.

33 K. Chiba, K. Yasuda and K. Nakamura (2001), "Numerical solution of fiber suspension flow through a parallel plate channel by coupling flow field with fiber orientation distribution", *Proc. 3rd Pacific Rim Conference on Rheol.*, Vancouver.
34 F. Chinesta and G. Chaidron (2001), "On the steady solution of linear advection problems in steady recirculating flows", *J. Non-Newtonian Fluid Mech.*, **98**, 65.

35 F. Chinesta, T. Mabrouki and A. Ramon (2002), "Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes", *5th Esaform Conference*.

36 S. Chono and M. Makino (1995), "Numerical simulation of fiber suspension flow between parallel plates", *Trans-B, Mach. Soc. Japan*, **61**, 3190.

37 R.M. Christensen (1990), "A critical evaluation for a class of micro-mecanics models", *J. Mech. Phys. Solids*, **38**, 379.

38 S.T. Chung and T.H. Kwon (1995), "Numerical simulation of fiber orientation in injection molding of short-fiber-reinforced thermoplastics", *Polym. Eng. Sci.*, **35**, 604.

39 T. Coupez, D. Daboussy and E. Bigot (1999), "Mesh adaptation for 3D injection numerical simulation", *15th Polymer Processing Society Conference*.

40 J.W. Daily and G. Bugliarello (1961), "Basic data for dilute fiber suspensions in uniform flow with shear", *TAPPI*, **44**, 497.

41 S.M. Dinh and R.C. Armstrong (1984), "A rheological equation of state for semiconcentrated fiber suspensions", *J. Rheol.*, **28**, 207.

42 M. Doi (1981), "Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropic and liquid crystalline phases", *J. Polym. Sci., Polym. Phys. Ed.*, **19**, 229.

43 A.A. Draad and M.A. Hulsen (1995), "Transition from laminar to turbulent flow for non-Newtonian fluids", in *Advances in Turbulence V*, R. Benzi (ed.), 105.

44 P.G. Drazin and W.H. Reid (1981), *Hydrodynamic Stability*, Cambridge University Press.

45 H.D. Ellis (1970), "Effects of shear treatment on drag reducing polymer solutions and fiber suspensions", *Nature*, **226**, 352.

46 J.L. Ericksen (1960a), "Anisotropic fluids", *Arch. Rat. Mech. Anal.*, **4**, 231.

47 J.L. Ericksen (1960b), "Transversely isotropic fluids", *Kolloidzeitschrift*, **173**, 117.

48 J.L. Ericksen (1966), "Instability in Couette flow of anisotropic fluids", *Quart. J. Mech. App. Math.*, **19**, 455.

49 K.A. Ericsson, S. Toll and J.A.E. Manson (1997), "The two-way interaction between anisotropic flow and fiber orientation in squeeze flow", *J. Rheol.*, **41**, 491.

50 J.D. Eshelby (1957), "The determination of the elastic field of an ellipsoidal inclusion and related problems", *Proc. Roy. Soc. London, A241*, 376.

51 L.G.R. Filipsson, J.H. Torgny Lagerstedt and F.H. Bark (1977), "A note on the analogous behaviour of turbulent jets of dilute polymer solutions and fibre suspensions", *J. Non-Newtonian Fluid Mech.*, **3**, 97.

52 G.E. Gadd (1965), "Friction reduction characteristics of fibrous and colloidal substances", *Nature*, **206**, 463.

53 P.G. de Gennes (1990), *Introduction to Polymer Dynamics*, Cambridge.

54 P. Gilormini and P. Vernusse (1992), "Tenseur d'Eshelby et problème d'inclusion dans le cas isotrope transverse incompressible", *C.R. Acad. Sci. Paris*, **314**, 257.

55 M. Grmela and P.J. Carreau (1987), "Conformation tensor rheological models", *J. Non-Newtonian Fluid Mech.*, **23**, 271.
56 V. Gupta, R. Sureshkumar, B. Khomami and J. Azaiez, "Linear stability of Taylor-Couette flow of semi-dilute non-Brownian fiber suspensions", submitted to *Phys. Fluids*.

57 M.F. Hibberd, M. Kwade and R. Scharf (1982), "Influence of drag reducing additives on the structure of turbulence in a mixing layer", *Rheol. Acta*, **21**, 582.

58 E.J. Hinch and L.G. Leal (1975), "Constitutive equations in suspension mechanics. Part 1", *J. Fluid. Mech.*, **71**, 481.

59 E.J. Hinch and L.G. Leal (1976), "Constitutive equations in suspension mechanics. Part 2: Approximate forms for a suspension of rigid particles affected by Brownian rotations", *J. Fluid Mech.*, **76**, 187.

60 E.J. Hinch (1977), "Mechanical models of dilute polymer solutions in strong flows", *Phys. Fluids*, **20**, S22.

61 J.W. Hoyt (1972), "The effect of additives on fluid friction", *ASME, J. Basic Eng.*, **94D**, 258.

62 G.B. Jeffery (1922), "The motion of ellipsoidal particles immersed in viscous fluid", *Proc. R. Soc., A102*, 161.

63 D.D. Joseph (1990), *Fluid Dynamics of Viscoelastic Liquids*, Springer Verlag.

64 D.D. Kale and A.B. Metzner (1976), "Turbulent drag reduction in dilute fiber suspensions: Mechanistic considerations", *AICHE J.*, **22**, 669.

65 D.L. Koch (1995), "A model for orientational diffusion in fibre suspensions", *Phys. Fluids*, **7**, 2086.

66 I.M. Krieger and T.J. Dougherty (1959), "A mechanism for non Newtonian flow in suspension of rigid spheres", *Trans. Soc. Rheol.*, **3**, 137.

67 M.T. Landahl (1973), "Drag reduction by polymer addition", In *Proceedings of Thirteenth International Congress of Theoretical and Applied Mechanics*, E. Becker and G.K. Mikhailov (Eds.), 177, Springer Verlag.

68 W.K. Lee, R.C. Vaseleski and A.B. Metzner (1974), "Turbulent drag reduction in polymeric solutions containing suspended fibers", *AICHE J.*, **20**, 128.

69 F.M. Leslie (1964), "The stability of Couette flow of certain anisotropic fluids", *Proc. Camb. Phil. Soc.*, **60**, 949.

70 F.M. Leslie (1966), "Some constitutive equations for anisotropic fluids", *Quart. J. Mech. Appl. Math.*, **XIX**, 357.

71 F.M. Leslie (1968), "Some constitutive equations for liquid crystals", *Arch. Rat. Mech. Anal.*, **28**, 265.

72 G.G. Lipscomb (1986), Ph.D dissertation, University of California, Berkeley.

73 G.G. Lipscomb, M.M. Denn, D.H. Hur and D.V. Boger (1988), "The flow of fiber suspensions in complex geometries", *J. Non-Newtonian Fluid Mech.*, **26**, 297.

74 J.L. Lumley (1969), "Drag reduction by additives", *Ann. Rev. Fluid Mech.*, **1**, 367.

75 W.D. McComb and K.T. Chan (1979), "Drag reduction in fiber suspensions: transitional behaviour due to fibre degradation", *Nature*, **280**, 45.

76 F. Meslin (1997), *Propietés rhéologiques des composites fibres courtes à l'état fondu*, PhD Thesis, ENS. Cachan.

77 F. Meslin, A. Poitou and F. Chinesta (1998), "Identification du comportement rhéologique de suspensions de fibres courtes", *Les Cahiers de Rhéologie*, **16**, 209-216.
78 F. Meslin, A. Poitou and F. Chinesta (1999), "A new rheometer for molten fiber composites", *15th Polymer Processing Society Conference*.

79 W. Mih and J. Parker (1967), "Velocity profile measurement and a phenomenological description of turbulent fiber suspension pipe flow", *TAPPI*, **50**, 237.

80 I. Monton (2002), *Numerical simulation of hort fiber composites forming processes*, PhD Thesis, University of Valencia (Spain).

81 T. Mura (1982), *Micromechanics of Defects in Solids*, Martinus Nijhoff.

82 B. Nsom (1994a), "Transition from circular Couette flow to Taylor vortex flow in dilute and semi-concentrated suspensions of stiff fibers", *J. Phys.*, **II 4**, 9.

83 B. Nsom (1994b), "Computation of drag reduction in fiber suspensions", *Fluid Dyn. Res.*, **14**, 275.

84 B. Nsom (1996), "Stability of fiber suspension flow in curved channel", *J. Phys.*, **II 6**, 1483.

85 T.C. Papanastasiou and A.N. Alexandrou (1987), "Isothermal extrusion of non-dilute fiber suspensions", *J. Non-Newtonian Fluid Mech.*, **25**, 313.

86 N. Phan-Thien and A.L. Graham (1990), "The squeezing flow of a model suspension fluid", *Rheol. Acta*, **29**, 433.

87 N. Phan-Thien and A.L. Graham (1991), "A new constitutive model for fibre suspensions: flow past a sphere", *Rheol. Acta*, **30**, 44.

88 E. Pichelin and T. Coupez (1998), "Finite element solution of the 3D mold filling for viscous incompressible fluid", *Comput. Meth. Appl. Mech. Engng.*, **163**, 359.

89 V.N. Pilipenko, N.M. Kalinichenko and A.S. Lemak (1981), "Stability of the flow of a fibre suspension in the gap between coaxial cylinders", *Sov. Phys. Dokl.*, **26**, 646.

90 F.T. Pinho and J.H. Whitelaw (1990), "Flow of non-Newtonian fluids in a pipe", *J. Non-Newtonian Fluid Mech.*, **34**, 129.

91 R.J. Pirih and W.M. Swanson (1972), "Drag reduction and turbulence modification in rigid particle suspensions", *Can. J. Chem. Eng.*, **50**, 228.

92 O. Pironneau (1989), *Finite methods for fluids*, Wiley.

93 A. Poitou, F. Chinesta and R. Torres (2000), "Numerical simulation of the steady recirculating flows of fiber suspensions", *J. Non-Newtonian Fluid Mech.*, **90**, 65.

94 M. Rahnama, D.L. Koch and E.S.G. Shaqfeh (1995), "The effect of hydrodynamic interactions on the orientation distribution in a fibre suspension subject to simple shear flow", *Phys. Fluids*, **7**, 487.

95 S. Ranganathan and S.G. Advani (1993), "A simultaneous solution for flow and fiber orientation in axisymmetric diverging radial flow", *J. Non-Newtonian Fluid Mech.*, **47**, 107.

96 J. Rosenberg, M.M. Denn and R. Keunings (1990), "Simulation of non-recirculating flows of dilute fiber suspensions", *J. Non-Newtonian Fluid Mech.*, **37**, 317.

97 S. Sasaki (1991a), "Drag reduction effect of rod-like polymer solutions: I. Influences of polymer cone concentration and rigidity of skeletal back bone", *J. Phys. Soc. Japan*, **60**, 868.

98 S. Sasaki (1991b), "Drag reduction effect of rod-like polymer solutions: II. Comparison between microgel and linear type polyions", *J. Phys. Soc. Japan*, **60**, 2613.

99 E.S.G. Shaqfeh and G.H. Frederickson (1990), "The hydrodynamic stress in a suspension of rods", *Phys. Fluids A*, **2**, 7.
100 E.S.G. Shaqfeh (1996), "Fully elastic instabilities in viscometric flows", *Ann. Rev. Fluid Mech.*, **28**, 129.

101 R.S. Sharma, V. Seshadri and R.C. Malhotra (1978), "Turbulent drag reduction by injection of fibers", *J. Rheol.*, **22**, 643.

102 M. Sussman, P. Smereka and S. Osher (1994), "A level set approach to computing solutions to incompressible two-phase flows", *J. Comput. Phys*, **114**.

103 W.G. Tiederman (1990), "The effect of dilute polymer solution on viscous drag and turbulent structure", In *Structure of Turbulence and Drag Reduction IUTAM Symp.*, A. Gyr (Ed.), Springer Verlag.

104 C.L. Tucker III (1991), "Flow regime for fiber suspensions in narrow gaps", *J. Non-Newt. Fluid. Mech.*, **39**, 239.

105 R.C. Vaseleski and A.B. Metzner (1974), "Drag reduction in the turbulent flow of fiber suspensions", *AICHE J.*, **20**, 301.

106 B.E. VerWeyst, C.L. Tucker III, P.H. Foss and J.F. O'Gara (1999), "Fiber orientation in 3-D injection molded features: prediction and experiment", *15th Polymer Processing Society Conference*.

107 P.S. Virk and D.L. Wagger (1990), "Aspects of mechanisms of type B drag reduction" In *Structure of Turbulence and Drag Reduction, IUTAM Symp.*, A. Gyr (Ed.), Springer Verlag.

108 Y. Yang, J.N. Chung, T.R. Troutt and C.T. Crowe (1990), "The influence of particles on the spatial stability of two-phase mixing layers", *Phys. Fluids*, **2**, 1839.