State-of-the-Art on numerical simulation of fiber-reinforced thermoplastic forming processes

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


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. 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. 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. surements (Meslin, 1997) 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. 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). 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   (Lipscomb et al. (1988)) 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. 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.

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 andLeal (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, axisymmetrical 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.

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 nondilute 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.
•  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 forthorder 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 equation 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 equation 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 orientation 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 considerable 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 tensors 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.

MECHANICAL MODELING: CONSTITUTIVE EQUATION
Two main steps allow to derive macroscopic constitutive equations for anisotropic suspensions: (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 Ω be a representative volume of the macroscopic scale, which contains many particles located in Ω i . Let τ , u and D denote respectively the microscopic stress tensor, velocity vector and strain rate tensor and n be a normal vector to the particle boundary ( Figure 9). For a given fiber orientation p (unit vector aligned along the principal axis of the fiber) in the suspension, for a Newtonian ambient fluid of viscosity η and if inertia terms can be neglected, the corresponding macroscopic variables T, v, and D are related by In the above expressions, p is the volume averaged pressure, I is the unit tensor and p represents the contribution of the particle to the averaged extra stress tensor.
Statistical average. The fiber orientation p is given in an average sense by its probability distribution ψ( p) so that the macroscopic statistical averaged stress tensor (i.e. the stress tensor used in engineering computations) is For most engineering applications, the exact calculation of int ψ( p) at each point is neither necessary nor possible. This orientations information is correct ly approximated with the use of the second order and fourth order orientation tensors, defined as follows: where ⊗ 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 λ denotes the aspect ratio of the particle, and if Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω and D are respectively the volume averaged vorticity and strain rate tensors, then The evolution equation of the probability distribution is given by 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 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 a 4 . The evolution equation for the fourth order tensor involves the 6 th order tensor, and in general, the evolution equation for any order moment of p involves the next higher-order moment. Therefore, to obtain a closed set of equations, one needs to determine a 4 as a function of 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 a 4 is defined as: (2.9)

Explicite Calculations for p
The calculation of p as a function of the macroscopic velocity field 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 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 (Batcherlor (1970-a)) is deduced, which has been modified by Shaqfeh and Fredrickson (1990). All of these models can be written in the following form In this expression, η 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 η eq ). In this case, if ϕ denotes the volume fraction of spheres, and if there are no hydrodynamic interactions (dilute sus-pension), the dilute approximation leads to Einstein's formula p = 5ηϕD( v) ⇒ η eq = η(1 + 2.5ϕ) (2.11) To account for hydrodynamic interactions, a natural way consists of calculating p by assuming that, instead of being immersed in the ambient fluid of viscosity η, the consequences of hydrodynamic interactions are summarized by placing the sphere in an ambient fluid of viscosity η eq (self-consistent model) (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 However the self-consistent model suffers from strong limitations because the equivalent viscosity increases to infinity when ϕ = 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 η eq (ϕ). We then add a quantity dϕ 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 ϕ, so that the volumetric rate induced by dϕ is dϕ/(1 − ϕ) 5 (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.

Algebraic preliminary result
A linear operator on the subspace of deviatoric (traceless) tensors, which exhibit a transverse isotropic symmetry of axis p depends on three parameters. A basis of such operators can be defined as This basis has two interesting properties. It exhibits the direction of transverse isotropy p and it allows easy products because it can be checked that:

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( v) = M is the constitutive relation of the fluid:M = 1 2η in the isotropic case, or more generally, M = 1 in the transverse isotropic case, then, the solution to generalized Jeffery's probem reads: 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 p is calculated by neglecting the particles interaction. Eq. (2.17) is thus solved with α 1 = α 2 = α 3 = 2η. 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 α 1 , α 2 , α 3 . With the same argument as the one that yielded to Eq. (2.12), it is possible to derive the following implicit equations: 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) : with initial conditions: α i = 2η for ϕ = 0.

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).
η 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 η 1 N p is equal to the parameter η 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 λ of 350, neither Shaqfeh and Fredrickson's nor Dinh and Armstrong's model is correct because, for very low values of ,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 account 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 selfconsistent 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.
where v is the velocity field.
• The constitutive equation for fibers with a large aspect ratio, and where the Brownian effects are neglected, is where p is the pressure field, I the unit tensor, D the strain rate tensor, a 4 the fourth order orientation tensor, η 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. (a 4 : D) ij = a 4 ijkl 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 where p is the unit vector aligned on the fiber axis direction, ψ( p) the orientation probability function and the symbol ⊗ denotes the tensorial product, i.e.
The fourth order orientation tensor is defined in a similar way by • Conservation balance of the fiber orientation distribution If we denote by p 0 the orientation of a fiber located at time t = 0 in the point x 0 , then the orientation of this fiber at time t, when it is located in the point x = x( x 0 , t) , will be referred to as p. The fiber orientation distribution at time t can be deduced from the initial one from the following conservation balance where ψ 0 ( p) denotes the initial fiber orientation distribution and | dp 0 | an infinitesimal solid angle defined around p 0 . To clarify the notation we consider a planar fiber distribution, where the fiber orientation can be defined by the angle ϕ with respect to a reference axis. Thus, the previous balance can be rewritten as where dϕ denotes the infinitesimal angle defined around the angle ϕ In general the initial fiber orientation distribution is known, and as we are going to show, the final orientation p can be computed from the initial one p 0 , i.e. we can compute the orientation angle ϕ of a fiber located at time t in the point x, from its initial orientation angle ϕ 0 . Thus, we can also compute the value of dϕ from dϕ 0 .
The orientation probability associated with the angle ϕ can be finally computed from • Fiber orientation transport In order to define the relation between the final and initial orientations ϕ = ϕ(ϕ 0 , t) we start by introducing Jeffery's equation which governs the orientation evolution of an ellipsoidal fiber immersed in a Newtonian fluid where k is a scalar parameter depending on the fiber aspect ratio λ (ratio between the fiber length and diameter) and Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω 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 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 The previous equations define the flow model of a suspension of short fibers. The rheological parameters η 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 = 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.

Boundary conditions
We denote by Ω f t the domain occupied by the fluid at time t. The boundary of the fluid domain which will be referred to as ∂Ω f t can be divided in two parts. In the first part, denoted by ∂ 1 Ω f t , the velocity is assumed known, whereas in the second one, ∂ 2 Ω f t , the traction is imposed. Moreover, the boundary ∂ 1 Ω f t can be divided in other two parts: the inflow boundary Γ − (injection gate, which is considered fixed in time) and the wall boundaries Γ 0 t where the adherence condition impose a zero velocity. In injection simulations, the boundary ∂ 2 Ω f t 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:

Fluid Domain Evolution: Flow Front Tracking
As commented, the flow model is defined in the part of the whole domain (injection mold) Ω occupied by the fluid at each time, t, Ω f t . In order to update the fluid domain we are going to introduce the fluid presence function, also known as fluid fraction function I( x, t). This function takes a unit value in the fluid region and it is zero in the empty domain The evolution of this function is given by the following scalar and linear advection equation This equation is defined in the whole domain. The fluid presence function must verify a boundary condition on the mold inlet (inflow boundary) 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

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 Ω f t = Ω.

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 v = 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:

and where the expression of the stress tensor is given by the constitutive equation (3.3).
In the previous variational formulation H 1 (Ω) and L 2 (Ω) denote the usual Sobolev and Lebesgue functional spaces and H 1 0 (Ω) is the functional space of velocities in H 1 (Ω) vanishing on the domain boundary ∂Ω. The election of the functions f (I), α(I) and β(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 α(I) and β(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 α and β were considered in that work. To justify the required values of both parameters, we consider an element Ω e which just starts its filling process. In this case the fluid fraction in this element takes an intermediate value 0 < I e 1. If the parameters α and β 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 α and β are too small, the equation of motion as well as the fluid incompressibility will be over-imposed in the element Ω 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 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 ≥ I th the flow model will be enforced in the element Ω e . However, as we shall show later, a judicious choice of the parameters α and β 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 Ω e a generic element where a linear approximation of the velocity field is considered. Thus, the velocity in each point of Ω 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 ≥ 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 Ω e at time t. We will assume that the element located upstream from Ω e , Ω 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 Ω 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 Ω 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 α and β (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 α(I e ) = 1 I e = 0 ε α I e > 0 (3. 24) and where ε α and ε β are two constants small enough, the nodal velocity at node x e is dominated by the flow equations imposed in the element Ω e−1 (the contribution of the pseudo-behavior imposed in the element Ω e (0 < I e < I th ), and consequently affected by the constants ε α and ε β 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 pseudobehavior (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 that 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 φ whose zero value determines the interface position. The evolution of this function is given by the following linear advection equation where u is an arbitrary velocity field which must coincide with the material velocity v on the interface at time t, and which will denoted from now on by Γ i (t): 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, φ( x, t), as well as the flow kinematics v( x, t) are known, then the interface at this time is given by (3.28) The solution of Eq. (3.26) allows to compute the level function at time t + Δt, φ( x, t + Δt), from which we can update the interface position (3.29) In the above equation, Δ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 φ( 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 where the functionφ coincides with φ along the interface, i.e.
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 mold. 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.

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 Taking into account the fluid incompressibility, Div v = 0, Eq (3.32) can be rewritten in the form dI dt = ∂I ∂t + Div( vI) = 0 (3.33) Now, we can consider the integral conservation form related to the previous differential equation where ∂Ω e denotes the boundary of Ω e , and n is the outwards unit vector defined on the element boundary. The element boundary can be divided in two parts: ∂Ω e + through which the flow leaves the element Ω e (i.e. the part of the element boundary verifying the relation v n > 0); and the inflow element boundary ∂Ω e − , though which the fluid is coming to the element Ω e (verifying v n < 0). Thus, the mass balance results 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 Ω e is approximated taking into account the fluid fraction existing in the upstream element Ω e− , I e− . In the same way, the outflow will depend on the fluid fraction existing in the element Ω e . Thus, the previous equation can be written as where |Ω e | is the volume of the element Ω e , q e =  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 Ω 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 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.

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. (1998Meslin et al. ( ), (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).

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.

Governing equations
The motion of an incompressible material is governed by mass and momentum conservation laws: where v is the velocity vector, ρ is the density of fluid, p is the isotropic pressure, τ τ τ τ τ τ τ τ τ τ τ τ τ τ 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 τ τ τ τ τ τ τ τ τ τ τ τ τ τ = 2ηD + φμD : a 4 (4.3) where η is the viscosity of the suspending Newtonian fluid, D is the deformation rate tensor, φ is the volume fraction of fibers, μ 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, μ reduces to where λ is the fiber aspect ratio (length to diameter ratio). The computations of fiber suspension flows through a 4:1 circular contraction (contraction ratio β = 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 p = v/ 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, φμ/η = φλ 2 / ln(λ), was varied, which can be determined completely by the volume fraction and aspect ratio of the fibers, φ and λ, respectively. It is called the fiber parameter.

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 μm and length l = 3 mm resulting in an aspect ratio λ = 231, and VF2 of diameter d = 11 μm and length l = 1 mm leading to λ = 94. These fibers were suspended in 65 wt% corn syrup/water solution of viscosity was η = 0.08 Pa·s at 16.5 o 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 φμ/η = 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 φμ/η = 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)). Furthermore, Figure 16 shows the fiber parameter φμ/η 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 φμ/η, 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 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 Figure8 in Chiba et al. (1994-b)). 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 φμ/η = 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η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

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.

Governing equations
The evolution equation of the orientation unit vector 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 The parameter k is given by then the orientation state can be illustrated using an orientation ellipse shown in Figure  19(b). The eigenvalues and the eigenvectors of 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. 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 This equation involves the fourth-order orientation tensor a 4 , then a closure equation is required to approximate a 4 using a. We used a quadratic closure approximation given by a 4 = a ⊗ a (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 ∂/∂t in a Lagrangian sense, then the evolution equations of the components of a may be written as ∂a 11 ∂t = 2Ω 12 + a 12 + 2k(D 11 a 11 + D 12 a 12 − D 11 a 2 11 − 2D 12 a 11 a 12 − D 22 a 11 a 22 ) (4.10) ∂a 12 ∂t = Ω 12 (a 22 − a 11 ) + k(D 12 − 2D 11 a 11 a 12 − 4D 12 a 2 12 − 2D 22 a 12 a 22 ) (4.11) The second method is a statistical scheme: orientation of a large number of fibers (N = 3600) with aspect ratio λ = 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 β is an angle between the x-axis and the fiber axis which is measured counterclockwise.
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 λ = 5 and 50, and the initial orientation was isotropic.

Performance of a quadratic closure approximation
The evolution of the orientation distribution function and the orientation ellipses for λ = 5 along Ψ = −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 Ψ 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 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.

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

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 λ = 5 and 10000 (k ≈ 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 = ρU H/η.

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 μm in diameter, l = 3 mm in length (λ = 283) were suspended in 63 wt% corn syrup/water solution with volume fraction of φ = 2 × 10 −5 . The viscosity of the Newtonian solvent is η = 0.0463 Pa·s at 16.8 • C. In this case the flow fields were almost the same as those for the solvent because the fiber parameter φμ/η = 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.  Figures 23 and 24 show fiber orientation state along three streamlines in a recirculating flow at Re = 15 for λ = 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.

Fiber orientation in a recirculating flow
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. 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 aspectratio 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 ψ * = 1.08 streamline for λ = 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 ψ * = 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).

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

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.

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 = ρU H/η. 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. where β is an angle between the x-axis and the projection of the fiber in the flow plane (x − y plane), and θ 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

Velocity profile and fiber orientation
Comparison of centerline velocity profile between the planar orientation and the threedimensional 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 ∼ 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 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 threedimensional 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.

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 ≤ x * ≤ 2 in Figures 31 and 32, respectively. The non-dimensional stress is defined as T * ij = T ij /(ηU/H). In the fully developed Newtonian flow of the channel T * xy is varied linearly from 6 at the wall to 0 at the centerline in the width direction, while T * xx − T * yy becomes zero. In contrast, complex distribution of stress components for fiber suspensions are remarkably observed in the region near the inlet owing to the orientation distribution of fibers. Both T * xy and T * xx − T * yy increase significantly from the values for Newtonian flow in the inlet near the wall. Furthermore, the negative value of T * xx − T * yy becomes larger in the inlet region near the centerline: deceleration of a fluid element produces a large compression along the centerline near the inlet, while an extension takes place near the wall and T * xx − T * yy becomes large. The effect of the three-dimensional orientation of fibers on the stress field are as follows: the degree of isotropy for the three-dimensional orientation is stronger than that for the planar orientation case, then, the variations of the stresses for the three-dimensional orientation of fibers are less significant than those for the planar case, such as the maximum values of T * xy are 14.4 and 11.1 for the planar and three-dimensional orientation, respectively. Furthermore, The maximum and minimum values of T * xx − T * yy are 13.9, -3.44 and 10.2, -1.82, respectively.

Conclusions
We studied the development of two-dimensional fiber suspension flows through a rectangular channel by solving the coupled flow kinematics and fiber orientation distribution for the planar and three-dimensional orientation of fibers. A fully-developed flow with a parabolic velocity profile and isotropic fiber orientation were specified at the inlet of the channel. 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.

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

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 exiting 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 1977and 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.

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: 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 where τ is the shear stress, v the bulk flow velocity and ρ 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 andWhitelaw 1990, andTiederman 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 longitudinal 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).

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 ω 1 and ω 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 Ericksen (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 orientation 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