Numerical simulation of flow kinematics for suspensions containing continuously dispersed aspect ratio fibers

Fiber suspension flow and fiber orientation through a parallel-plate channel were numerically simulated for fiber suspensions including continuously dispersed aspect ratios from 10 to 50. In the simulations, both the fiber–fiber and fiber–wall interactions were not taken into account. A statistical scheme that proceeds by evaluating the orientation evolution of a large number of fibers from the solution of the Jeffery equation along the streamlines was confirmed to be a very useful and feasible method to accurately analyze the orientation distribution of fibers with continuously dispersed aspect ratios. For monodisperse suspensions with small-aspect-ratio fibers, flip-over or oscillation phenomenon of the orientation ellipsoid caused the wavy patterns of the velocity profile and the streamlines as well as the abrupt and complex variation of the shear stress and the normal stress difference near the channel wall as proven in one of our former works. On the other hand, continuous dispersions containing from small- to large-aspect-ratio fibers were able to induce smoother evolutions of the fiber orientation and the flow kinematics. In the processing of fiber composites, the length of suspended fibers is always continuously distributed because of fiber breakage during processing; thus, the smooth evolutions of the flow kinematics and the stress distribution can be attained.


Kunji Chiba Francisco Chinesta
Numerical simulation of flow kinematics for suspensions containing continuously dispersed aspect ratio fibers Abstract Fiber suspension flow and fiber orientation through a parallelplate channel were numerically simulated for fiber suspensions including continuously dispersed aspect ratios from 10 to 50. In the simulations, both the fiber-fiber and fiber-wall interactions were not taken into account. A statistical scheme that proceeds by evaluating the orientation evolution of a large number of fibers from the solution of the Jeffery equation along the streamlines was confirmed to be a very useful and feasible method to accurately analyze the orientation distribution of fibers with continuously dispersed aspect ratios. For monodisperse suspensions with smallaspect-ratio fibers, flip-over or oscillation phenomenon of the orientation ellipsoid caused the wavy patterns of the velocity profile and the streamlines as well as the abrupt and complex variation of the shear stress and the normal stress difference near the channel wall as proven in one of our former works. On the other hand, continuous dispersions containing from small-to large-aspect-ratio fibers were able to induce smoother evolutions of the fiber orientation and the flow kinematics. In the processing of fiber composites, the length of suspended fibers is always continuously distributed because of fiber breakage during processing; thus, the smooth evolutions of the flow kinematics and the stress distribution can be attained.

Introduction
Fiber suspension can often exhibit anisotropy due to flowinduced fiber alignment. The addition of short fibers to a Newtonian liquid can therefore drastically change the flow kinematics even at very low concentrations. Furthermore, the fiber orientation state strongly affects the mechanical properties of short-fiber composites, and the characterization of fiber orientation distribution becomes much more significant in current processing of high-quality composite materials. The importance of coupling flow kinematics and fiber orientation, therefore, has been recognized during the last two decades (Ausias et al. 1992(Ausias et al. , 1994Azaiez et al. 1997;Chiba et al. 1990;Chung and Kwon 2002;Fan et al. 1999;Ghosh et al. 1995;Lipscomb et al. 1988;Rosenberg et al. 1990;Zirnsak and Boger 1998), and we can also find extensive literature focusing on fiber suspension flow in Chiba and Chinesta (2005).
In a previous paper, Chiba and Chinesta (2002) simulated the 2-D suspension flows involving both quasi-infinite and small-aspect-ratio (r a =10) fibers in a Newtonian fluid through a parallel-plate channel by coupling flow kinematics with the 3-D fiber orientation distribution. The influence of the flip-over phenomenon on the flow kinematics and fiber orientations was noted. It is well known that fibers with small aspect ratio frequently and periodically flip over when fiber-fiber interaction can be neglected. Flip-over is a phenomenon related to the rapid rotation of a slender particle K. Chiba Faculty of Education, Shiga University, 2-5-1 Hiratsu, Otsu, Shiga 520-0862, Japan e-mail: kchiba@sue.shiga-u.ac.jp F. Chinesta LMSP: UMR CNRS-ENSAM-ESEM, 151 Boulevard de l'Hôpital, 75013 Paris, France e-mail: francisco.chinesta@paris.ensam.fr that occurs when it is aligned in the direction perpendicular to the flow. The just referred study focused on monodisperse fiber suspensions.
In our recent work (Chiba and Chinesta 2005), the development of flow kinematics and fiber orientation distribution through a parallel-plate channel was simulated for multidisperse suspensions with discretely dispersed aspect ratio fibers, and their predictions were compared with those obtained for mono-and bidisperse suspensions. For suspensions including small-aspect-ratio fibers, flip-over or oscillation phenomenon of the orientation ellipsoid (schematic illustration of fiber orientation distribution) caused the wavy patterns of the velocity profile and the streamlines as well as the rapid and complex variation of the shear stress and the normal stress difference near the channel wall. In the processing of fiber composites, however, the length of suspended fibers is always continuously distributed because of fiber breakage during processing (Shiraishi et al. 2004) except for compression molding. Thus, the effect of continuous distribution of fiber lengths on the flow kinematics should be examined.
In this work, coupled flow kinematics and 3-D fiber orientation distribution were simulated for multidisperse suspensions with continuously distributed aspect ratio fibers; then the mechanism of the flip-over or oscillation phenomenon of the orientation ellipsoid was examined in detail.
Fiber suspensions that are the focus of the present work involve rigid fibers in a Newtonian liquid. However, in spite of the significant impact of fiber-fiber and fiber-wall interactions on the suspension behavior, in the present paper these interactions were not considered (the introduction of diffusion terms that account for these effects is a work in progress). A parabolic velocity profile and an isotropic orientation distribution were prescribed on the inflow boundary and a natural flow-out condition on the exit of a two-parallel-plate channel. The development of the flow kinematics and 3-D fiber orientation is analyzed.

Numerical solution procedure for multidisperse suspension flow
The coupled flow kinematics and 3-D fiber orientation distribution model was simulated to study the development of 2-D flows (flow plane is the xy plane) of suspensions with continuously dispersed fibers in a Newtonian fluid through a two-parallel-plate channel. Detailed description of the governing equations, statistical scheme, and solution procedure can be found for multidisperse fiber suspension in Chiba and Chinesta (2005), which we will briefly summarize in the next paragraphs.

Governing equations
Mass and momentum conservation equations govern the motion of an incompressible material: where v is the velocity vector, ρ is the fluid density, 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 to account for the effect of the presence of fibers on the stress field. The extra stress can be written as where η is the viscosity of the suspending Newtonian fluid, D is the strain rate tensor, φ is the volume fraction of fibers, μ is a material constant, and a 4 denotes the fourth-order orientation tensor defined by where p is the unit vector defining the fiber axis and ψ(p) is the orientation distribution function (the orientations of fibers in a fluid element, which flows along a streamline, are distributed). For the aspect ratio of fiber (length to diameter ratio) r a ≥ 10, μ reduces to (Lipscomb et al. 1988) Statistical scheme The evolution of p in a Newtonian homogeneous flow can be described by the Jeffery equation (Jeffery 1922), where Ω is the vorticity tensor and the parameter k is given by The second-order orientation tensor a 2 (see Eq. 8) introduced by Advani and Tucker (1987) is a suitable and concise way of describing the orientation distribution of fibers. Schematic orientation state can be illustrated using an orientation ellipsoid (e.g., Fig. 2). The eigenvalues and the eigenvectors of a 2 give the lengths and directions of the three major axes of that ellipsoid. The eigenvalues indicate the intensity of fiber orientation along the corresponding directions (eigenvectors).
We used a statistical scheme in the computations of fiber orientation distribution, which proceeds by evaluating the 3-D orientations of a large number of fibers with specified initial angles by solving the Jeffery equation (6) along the streamlines using a fourth-order Runge-Kutta technique. The appropriate number of segmentations between 0°and 90°in the azimuthal direction was 30 or more. In the present simulations, however, the number of segmentations was set to 18 (the number of fibers equals 816 for isotropic orientation) for each aspect ratio fiber to reduce the required CPU time. We were able to accurately estimate the orientation distribution of fibers by 18 segmentations (Chiba and Chinesta 2005). Then, the components of a 2 and a 4 can be calculated along the streamlines according to Eqs. (9) and (10) in Chiba and Chinesta (2005).

Solution procedure
All governing equations were written in terms of nondimensional variables using the channel width H and the mean fluid velocity U as the reference parameters. Subsequently, the Reynolds number was defined by Re ¼ UH=: The computational domain was x*=0 to 10 in the flow direction and y*=0 to 0.5 in the width direction (the lower-half part of the channel). We used the finite difference method for the flow kinematics resolution. The details concerning the numerical procedure used for computing the extra stress can be found in Eqs. (11) and (12) in Chiba and Chinesta (2005).
The simulations were carried out for the five scenarios of Table 1. For the bidisperse suspension and scenarios 1 to 3 concerning multidisperse suspensions, the volume fraction of each fiber aspect ratio included in the suspension is the same. On the other hand, the volume fraction in the fourth scenario is distributed as shown in Fig. 1. The histogram shows the volume fraction distribution of glass fibers in the injection-molded specimen measured by Shiraishi et al. (2004), and its fitting curve was used in the simulation. Thus, the volume fraction increases and reaches the maximum value for r a ≈23, then gradually decreases as the aspect ratio is further increased. Furthermore, the Reynolds number and the total fiber parameter were set to Re=10 and = ¼ 10; respectively.
For multidisperse suspensions, the orientation ellipsoids are determined from the orientation distribution of all fibers, including different aspect ratio fibers. Of course, the number of fibers related to each aspect ratio included in the suspension depends on the volume fraction and the fiber aspect ratio; for example, for r a =10+50 (bidisperse suspension), the number of fibers with r a =50 results in one fifth of that for r a =10 fibers, because the volume fraction being the same, the length ratio is 5.

Evolution of fiber orientation in multidisperse suspension flow
In this section, the mechanism associated with the flip-over or oscillation phenomenon of the orientation ellipsoid for multidisperse suspensions will be discussed.  We can see clearly the evolutions of the projections of the orientation ellipsoids on both the xy plane (flow plane; Fig. 2a) and the xz plane (Fig. 2b) along the streamlines related to the coordinates y*=1/40, 2/40 and 3/40 at the inlet (x*=0) of the two-parallel-plate channel for the bidisperse suspension.
The orientation ellipsoids frequently and periodically flip over or oscillate along the streamlines for suspensions, including small-aspect-ratio fibers, and also the flip-over or a projections on the xy-plane b projections on the xz-plane Fig. 2 Evolution of the orientation ellipsoids for the r a =10+50 (bidisperse suspension) along the streamlines related to the coordinates y*=1/40, 2/40, and 3/40 at the channel inlet; projections on the xy plane (a) and xz plane (b) a scenario 1 b scenario 2 c scenario 3 Fig. 3 Evolution of the orientation ellipsoids for multidisperse suspensions along the streamlines related to the coordinates y*=1/40, 2/40, and 3/40 at the channel inlet (projections on the xy plane) Fig. 4 Evolution of the second-order orientation tensor components and the director angle along the streamline defined by y*=1/40 at the channel inlet for bi-and multidisperse suspensions oscillation phenomenon increases as the channel wall is approached because of higher shear rates. Furthermore, Fig. 2b clearly shows that the preferential direction of the orientation ellipsoid (eigenvector related to the highest eigenvalue, also now as director) always lies on the xy plane even when the flip-over or oscillation phenomenon occurs because the initial orientation of fibers was isotropic.
Note that the director almost completely aligns along the streamline as the number of aspect ratios included in the suspension increases by comparing the evolutions of the fiber orientation for the bidisperse case (Fig. 2) and scenarios 1, 2, and 3 involving multidisperse suspensions (Fig. 3). Figure 4 illustrates the variations of the second-order orientation tensor components and the director angle for all the computational conditions. The director angle abruptly changes from a negative value to a positive one, and two types of abrupt changes in the director angle can be found: (a) it continuously changes through 0 (occurrence of the oscillation phenomenon) and (b) it reaches ±90°(occurrence of flip-over at x*=4). On the other hand, a 11 becomes minimum, while a 22 , maximum at the positions of occurrence of flip-over or oscillation phenomenon: flipover phenomenon occurs when the minimum value of a 11 becomes smaller than the maximum value of a 22 , i.e., a 11 <1/3. In contrast, oscillation phenomenon occurs at the positions where the minimum value of a 11 is larger than the maximum value of a 22 . Furthermore, note that the amplitudes of the fluctuation of a 11 , a 22 , and the director angle decrease as the number of aspect ratios increases. Figure 5 shows the detailed evolution of the orientation ellipsoid within the region of flip-over or oscillation phenomenon for the r a =10 monodispersion (I), the bidispersion (II), and the four multidisperse scenarios (III to VI). Black-filled ellipses clearly illustrate flip-over for the mono-and bidispersions and oscillation phenomenon for the multidispersions.
In Figs. 2, 3, 4, 5, the orientation ellipsoids involve all fiber aspect ratios. The evolutions of the orientation ellipsoids related to each fiber aspect ratio, e.g., for scenario 1, are shown in Fig. 6 to better understand the mechanism of the flip-over and oscillation phenomenon in multidisperse suspensions. The region shown in Fig. 6 (approximately 3.12≤x*≤3.18) corresponds to the region Fig. 5 Detailed evolution of the orientation ellipsoids for mono-, bi-, and multidisperse suspensions: r a =10, monodispersion (I), 1.4≤x*≤1.6; r a =10+50, bidispersion (II), 3.9≤x*≤4.1; scenario 1 (III), 3.05≤x*≤3.25; scenario 2 (IV), 4.8≤x*≤5.0; scenario 3 (V), 4.8≤x*≤5.0; and scenario 4 (VI), 4.8≤x*≤5.0 Fig. 6 Oscillation phenomenon of the orientation ellipsoid for scenario 1 (III, black-filled ellipses) and the evolution of the orientation ellipsoid for each fiber aspect ratio (r a =10, 20, 30, 40, and 50) illustrated by black-filled ellipses in Fig. 5, which is expressed by the symbol (III). The assembly of fibers with finite aspect ratio can align, flip over, or oscillate. In this specified region, however, fibers almost completely align in one direction for both r a =30 and 50, the orientation ellipsoid oscillates for r a =20 and 40 (gray-filled ellipse), and flip-over for r a =10 (gray-filled ellipse). The superposition of these evolutions of the orientation ellipsoid for each fiber aspect ratio included in the suspension results in the oscillation phenomenon when all the fibers involved in the suspension are considered to define the orientation ellipsoids (black-filled ellipses in Fig. 6).

Velocity profile and stress distribution in multidisperse suspension flow
We observe in Fig. 7a the wavy pattern of the centerline velocity profile for the r a =10 monodisperse suspension in contrast to the smooth profile for the quasi-infinite aspect ratio monodispersion. Furthermore, the centerline velocity profile gradually becomes smoother as the number of aspect ratios included in the suspension is further increased, as can be appreciated in the results obtained for the bidispersion and the four multidisperse scenarios. Under all computational conditions, the centerline velocity profile related to the fiber suspensions changes significantly with respect to the Newtonian flow near the 10.9 to 12.7 9.1 to 10.9 7.3 to 9.1 5.5 to 7.3 3.6 to 5.5 1.8 to 3.6 0.0 to 1.8 g quasi infinite aspect-ratio mono-dispersion Fig. 8 (continued) channel inlet. 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. The flow kinematics can reach a fully developed regime only when the number of aspect ratios included in the suspension increases as for the quasi-infinite aspect ratio monodispersion. The velocity profile of the x component in the width direction can be seen in Fig. 7b for scenario 4. It can be found that the velocity profile for suspension with smallaspect-ratio fibers becomes wavy except the region near the inlet of the channel even at the exit (see Fig. 6b in Chiba and Chinesta 2002 for the r a =10 monodisperse suspension), while the velocity profile becomes smoother as the number of aspect ratios included in the suspension increases.
Furthermore, Fig. 8a-f shows that the abrupt and complex distribution of the principal stress difference near the channel wall, which is observed for the r a =10 monodisperse suspension (Fig. 8a), is reduced as the number of aspect ratios increases. For scenario 4 (Fig. 8f), the distribution becomes very smooth and is very close to the one related to the quasi-infinite aspect ratio monodisperse case (Fig. 8g) because the intensity of flip-over or oscillation phenomenon can be suppressed and the orientation ellipsoid exhibits better alignment to the streamline as explained in the previous section.

Conclusion
The effects of continuously dispersed aspect ratio fibers on the flow kinematics and the fiber orientation were nu-merically analyzed in multidisperse fiber suspension flows, and the following conclusions were obtained: 1. Occurrence of flip-over or oscillation phenomenon of the orientation ellipsoid depends on the number of fibers of each aspect ratio as well as on the flow kinematics. 2. As the number of aspect ratios included in the suspension increases, the fiber assembly almost completely aligns along the streamline; thus, the velocity profile and the stress distribution become smooth, i.e., the wavy patterns of the centerline velocity profile and streamlines as well as the localized and complex variation of stresses near the channel wall are removed. Consequently, we can conclude that the mechanical and physical properties of composites should be improved by increasing number of aspect ratios included in the suspension as well as the volume fraction related to large fiber aspect ratios.
The suspension regime considered in the present simulations corresponds to a semiconcentrated or concentrated regime; thus, the Jeffery equation is no longer an appropriate way to describe the fiber orientation because it does not take account of the fiber-fiber interaction. This effect can be introduced in the Jeffery formalism from an appropriate Brownian term. This question will be addressed in future work. Furthermore, a still open problem that is encountered in the processing of composite materials is the modeling of the role of the interaction between fibers and polymer molecules on the evolution of the fiber orientation, i.e., an extension of the Jeffery equation to viscoelastic media.