Numerical Simulation of Flow Kinematics for
Suspensions Containing Continuously-Dispersed
Aspect-Ratio Fibers

Kunji Chiba, Francisco Chinesta

► To cite this version:

Kunji Chiba, Francisco Chinesta. Numerical Simulation of Flow Kinematics for Suspensions Containing Continuously-Dispersed Aspect-Ratio Fibers. Rheologica Acta, 2006, 45 (5), pp.547-554. 10.1007/s00397-005-0064-0 . hal-01007129

HAL Id: hal-01007129

https://hal.science/hal-01007129

Submitted on 5 Mar 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical simulation of flow kinematics for suspensions containing continuously dispersed aspect ratio fibers

Kunji Chiba
Francisco Chinesta

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

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

**Keywords** Multidisperse fiber suspension · Continuously dispersed aspect ratio fibers · Statistical scheme · Flip-over · Oscillation phenomenon

## Introduction

Fiber suspension can often exhibit anisotropy due to flow-induced 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, 1994; Azaiez 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
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:

$$Div \mathbf{v} = 0, \qquad (1)$$

$$\rho dv/dt = -Grad p + Div \tau, \qquad (2)$$

where $v$ is the velocity vector, $\rho$ is the fluid density, $p$ is the isotropic pressure, $\tau$ is the extra stress tensor, and $d/dt$ denotes the material derivative. The continuum theory for dilute suspensions of rigid, high-aspect-ratio particles developed by Lipscomb et al. (1988) is used to account for the effect of the presence of fibers on the stress field. The extra stress can be written as

$$\tau = 2\eta \mathbf{D} + \phi \mu \mathbf{D} : \mathbf{a}_4, \qquad (3)$$

where $\eta$ is the viscosity of the suspending Newtonian fluid, $\mathbf{D}$ is the strain rate tensor, $\phi$ is the volume fraction of fibers, $\mu$ is a material constant, and $\mathbf{a}_4$ denotes the fourth-order orientation tensor defined by

$$\mathbf{a}_4 = \oint \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}) d\mathbf{p}, \qquad (4)$$

where $\mathbf{p}$ is the unit vector defining the fiber axis and $\psi(\mathbf{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 \ge 10$, $\mu$ reduces to (Lipscomb et al. 1988)

$$\mu = \frac{\eta r_a^2}{\ln(r_a)}. \qquad (5)$$

### Statistical scheme

The evolution of $\mathbf{p}$ in a Newtonian homogeneous flow can be described by the Jeffery equation (Jeffery 1922),

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

where $\Omega$ is the vorticity tensor and the parameter $k$ is given by

$$k = \frac{(r_a^2 - 1)}{(r_a^2 + 1)}. \qquad (7)$$
The second-order orientation tensor $\mathbf{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 $\mathbf{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).

$$ \mathbf{a}_2 = \oint \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}) d\mathbf{p}. \qquad (8) $$

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 $\mathbf{a}_2$ and $\mathbf{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 non-dimensional variables using the channel width $H$ and the mean fluid velocity $U$ as the reference parameters. Subsequently, the Reynolds number was defined by $Re = \rho UH/\eta$. 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

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Volume fraction distribution of glass fibers in the injection-molded specimen; diameter of fiber $d=13$ µm, measured by Shiraishi et al. (2004)</figcaption></figure>

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 \approx 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 $\phi\mu/\eta = 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.

<table><caption>**Table 1** Computational conditions for bi- and multidisperse fiber suspensions</caption><thead><tr><th>Scenario</th><th>Aspect ratios of fibers ($r_a$)</th><th>Volume fraction for each fiber aspect ratio</th></tr></thead><tbody><tr><td>Bidispersion</td><td>10 and 50</td><td>The same</td></tr><tr><td>Scenario 1</td><td>10, 20, 30, 40, and 50</td><td>The same</td></tr><tr><td>Scenario 2</td><td>10, 15, 20, ..., 40, 45, and 50</td><td>The same</td></tr><tr><td>Scenario 3</td><td>10, 12, 14, ..., 46, 48, and 50</td><td>The same</td></tr><tr><td>Scenario 4</td><td>10, 12, 14, ..., 46, 48, and 50</td><td>Distributed as shown in Fig. 1</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**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)</figcaption></figure>

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.

<figure><img src="image_4.png" /><figcaption>**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)</figcaption></figure>

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

<figure><img src="image_5.png" /><figcaption>**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</figcaption></figure>
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: flip-over 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

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Detailed evolution of the orientation ellipsoids for mono-, bi-, and multidisperse suspensions: $r_a=10$, monodispersion (I), $1.4 \le x^* \le 1.6$; $r_a=10+50$, bidispersion (II), $3.9 \le x^* \le 4.1$; scenario 1 (III), $3.05 \le x^* \le 3.25$; scenario 2 (IV), $4.8 \le x^* \le 5.0$; scenario 3 (V), $4.8 \le x^* \le 5.0$; and scenario 4 (VI), $4.8 \le x^* \le 5.0$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**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, \text{ and } 50$)</figcaption></figure>

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 \le x^* \le 3.18$) corresponds to the region

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

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Velocity profiles in a developing flow of fiber suspensions through a two-parallel-plate channel; centerline velocity profiles for mono-, bi-, and multidisperse suspensions (a); velocity profile in the width direction for scenario 4 (b)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Principal stress difference in mono-, bi-, and multi-disperse suspension flows</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>e scenario 3 multi-dispersion</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>f scenario 4 multi-dispersion</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>g quasi infinite aspect-ratio mono-dispersion</figcaption></figure>

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 bi-dispersion 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
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 small-aspect-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 numerically 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.

## References

Advani SG, Tucker CL III (1987) The use of tensors to describe and predict fiber orientation in short fiber composites. J Rheol 31:751–784

Ausias G, Agassant JF, Vincent M, Lafleur PG, Lavoie PA, Carreau PJ (1992) Rheology of short glass fiber reinforced polypropylene. J Rheol 36:525–542

Ausias G, Agassant JF, Vincent M (1994) Flow and fiber orientation calculations in reinforced thermoplastic extruded tubes. Int Polym Process 9:51–59

Azaiez J, Guenette R, Ait-Kadi A (1997) Investigation of the abrupt contraction flow of fiber suspensions in polymeric fluids. J Non-Newton Fluid Mech 73:289–316

Chiba K, Chinesta F (2002) Simulation numérique et observation de l'écoulement d'une suspension de fibres dans une filière à géometrie complexe. Rev Compos Mater Av 12:499–513

Chiba K, Chinesta F (2005) Numerical simulation of flow kinematics and fiber orientation for multi-disperse suspension. Rheol Acta 45:1–13

Chiba K, Nakamura K, Boger DV (1990) A numerical solution for the flow of dilute fiber suspensions through an axisymmetric contraction. J Non-Newton Fluid Mech 35:1–14

Chung DH, Kwon TH (2002) Numerical studies of fiber suspensions in an axisymmetric radial diverging flow: the effects of modeling and numerical assumptions. J Non-Newton Fluid Mech 107:67–96

Fan X-J, Phan-Thien N, Zheng R (1999) Simulation of fiber suspension flows by the Brownian configuration field method. J Non-Newton Fluid Mech 84:257–274

Ghosh T, Grmela M, Carreau PJ (1995) Rheology of short fiber filled thermoplastics. Polym Compos 16:144–153

Jeffery GB (1922) The motion of ellipsoidal particles immersed in viscous fluid. Proc R Soc Lond A 102:161–179

Lipscomb GG, Denn MM, Hur DH, Boger DV (1988) The flow of fiber suspensions in complex geometries. J Non-Newton Fluid Mech 26:297–325

Rosenberg J, Denn MM, Keunings R (1990) Simulation of non-recirculating flows of dilute fiber suspensions. J Non-Newton Fluid Mech 37:317–345

Shiraishi Y, Narazaki N, Kikutani T (2004) Numerical predictions of 3-D filler orientation and elastic compliance for reinforced plastics containing fillers with aspect ratio distribution. Polym Compos 25:194–213

Zirnsak MA, Boger DV (1998) Axisymmetric entry flow of semi-dilute xanthan gum solutions: prediction and experiment. J Non-Newton Fluid Mech 79:105–136