# Modeling fiber interactions in semiconcentrated fiber suspensions

J. Férec, G. Ausias, M.C. Heuzey, P.J. Carreau

► To cite this version:

J. Férec, G. Ausias, M.C. Heuzey, P.J. Carreau. Modeling fiber interactions in semiconcentrated fiber suspensions. Journal of Rheology, 2009, 53, pp.49-72. 10.1122/1.3000732. hal-01552472

HAL Id: hal-01552472

https://hal.science/hal-01552472

Submitted on 3 Jul 2017

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

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

Public Domain
# Modeling fiber interactions in semiconcentrated fiber suspensions

J. Férec

Center for Applied Research on Polymers and Composites (CREPEC),
Chemical Engineering Department, École Polytechnique de Montréal,
PO Box 6079, Stn Centre-Ville, Montreal, Quebec H3C3A7, Canada
and Laboratoire d'Ingénierie des MATériaux de Bretagne (LIMATB),
Université Européenne de Bretagne, rue de St. Maudé, 56325 Lorient, France

G. Ausias

Laboratoire d'Ingénierie des MATériaux de Bretagne (LIMATB),
Université Européenne de Bretagne, rue de St. Maudé, 56325 Lorient, France

M. C. Heuzey and P. J. Carreau

Center for Applied Research on Polymers and Composites (CREPEC),
Chemical Engineering Department, École Polytechnique de Montréal,
PO Box 6079, Stn Centre-Ville, Montreal, Quebec H3C3A7, Canada

A set of rheological equations is developed for semiconcentrated suspensions of rigid fibers in a Newtonian fluid taking into account hydrodynamic and fiber-fiber interactions. The force generated by the fiber interactions is modeled using a linear hydrodynamic friction coefficient proportional to the relative velocity at the contact point, and weighted by the probability for contacts to occur. The equation of evolution of the second-order orientation tensor, containing advection and diffusion terms due to fiber interactions, is derived to predict fiber orientation under flow. The well known fourth-order orientation tensor, related to the hydrodynamic contribution, and a newly proposed fourth-order interaction tensor are used to evaluate the total stress in the composite. A linear and a quadratic closure approximation are proposed to describe the fourth-order interaction tensor. Results are presented using the quadratic form, which is found to be more accurate than the linear one. The model is shown to describe well simple shear data of suspensions of glass fibers in a Newtonian polybutene. Moreover, fiber orientation and the average number of contacts per fiber are predicted. The newly proposed interaction coefficient varies with fiber orientation, which appears to be realistic.

## I. INTRODUCTION

Thermoplastic polymer composites containing short fibers can be extruded or molded into the desired shape without any significant modification of conventional processing equipment. Owing to flow, the fiber orientation changes dynamically during processing and the final orientation pattern is retained in the part following solidification. From this
orientation state, it is possible to determine the mechanical and other physical properties of the final product. However, this is not a simple task since fiber motion and rotation are affected by fluid flow and vice versa.

Many investigations have been carried out to elucidate the rheological behavior of suspensions of elongated particles. Some pioneering theoretical developments have been presented by Jeffery (1922), who developed a model for a single dilute spheroid suspended in a Newtonian fluid. He calculated the stress around a particle in the absence of external torque and derived the evolution equation of the unit vector p directed along the particle principal axis.

Generally, fiber suspensions are concentrated systems and fibers cannot be considered as individual particles. Therefore, a reference cell that contains a large number of fibers needs to be defined [Batchelor (1971)]. The orientation state of this fiber population in the reference cell can be represented by a probability distribution function, from which Advani and Tucker (1987) have defined orientation tensors, second- and fourth-order tensors, $a_2$ and $a_4$, respectively. These describe more easily the fiber orientation state than the probability distribution function, without significant losses.

Fiber suspensions can be characterized according to the volumetric fraction of solid particles in the fluid. Depending on the fiber volume fraction $\phi$, some quantities such as hydrodynamic or fiber-fiber interaction forces dominate the system. Typically, the number of fibers per unit volume $n$ is divided into three distinct regimes [Doi and Edwards (1978)]. In the dilute regime defined by $n \ll 1/L^3$, where $L$ is the length of the particle, the fibers are allowed to move without interference from other particles. By increasing the number of fibers, $1/L^3 < n < 1/DL^2$ ($D$=particle diameter), the semidilute domain is reached and non-negligible hydrodynamic interactions between fibers are encountered. Furthermore, some fiber contacts are possible. Finally, the concentrated regime is characterized by the numerous contacts of a fiber with its neighbors, and in this state, $n > 1/DL^2$.

From experimental observations in concentrated systems the primary fiber orientation is that of the flow direction, but some deviations appear in injection parts [Bay (1991)] and extruded tubes [Ausias et al. (1994)]. Furthermore, Petrich et al. (2000) observed that their measured viscosity and first normal stress differences were much larger than the predictions of the hydrodynamic theories and attributed these enhancements to the presence of fiber-fiber contacts. The modeling of these effects represents a challenging task. Folgar and Tucker (1984) and then Advani and Tucker (1987) included into the Jeffery equation an empirical diffusional term that accounts for particle interactions. This phenomenological term is proportional to the effective deformation rate $|\dot{\gamma}| = \sqrt{\dot{\gamma}:\dot{\gamma}}/2$ and to a constant diffusion coefficient $C_I$, as shown below in the time evolution equation for the second-order orientation tensor [Advani and Tucker (1987)],

$$ \frac{D a_2}{Dt} = -\frac{1}{2}(\omega \cdot a_2 - a_2 \cdot \omega) + \frac{1}{2}\lambda(\dot{\gamma} \cdot a_2 + a_2 \cdot \dot{\gamma} - 2a_4:\dot{\gamma}) + 2C_I|\dot{\gamma}|(\delta - 3a_2), \quad (1) $$

where $\omega$ and $\dot{\gamma}$ are the vorticity and deformation rate tensors, respectively. $\delta$ is the identity tensor and the shape factor $\lambda = (r^2 - 1)/(r^2 + 1)$ is a parameter related to the particle aspect ratio $r=L/D$. Several authors have tried to express the diffusion coefficient $C_I$ as functions of the fiber aspect ratio and volume fraction. Based on experimental data, Bay (1991) proposed the following exponential function:

$$ C_I = 0.0184 \exp(-0.7148r\phi). \quad (2) $$
Phan-Thien et al. (2002) proposed the following expression based on calculations of the diffusion coefficient using a direct simulation of the fiber dynamics:

$$C_I = 0.03[1 - \exp(-0.224r\phi)]. \qquad (3)$$

Finally, Ranganathan and Advani (1991) assumed a simple relationship between the diffusion coefficient and the average interfiber spacing $a_c$ for semiconcentrated suspensions, with $C_I=KL/a_c$, where K is a constant.

In terms of constitutive equations, Hand (1961) and then Giesekus (1962) derived expressions for dilute ellipsoid suspensions in a Newtonian fluid of viscosity $\eta_0$. Lipscomb et al. (1988) have extended the previous work to propose a more general governing equation with corresponding material coefficients (A, B, C, and F) for dilute ellipsoid suspensions, represented by

$$\sigma = -P\delta + \eta_0\dot{\gamma} + \eta_0\phi[Aa_4:\dot{\gamma} + B(a_2 \cdot \dot{\gamma} + \dot{\gamma} \cdot a_2) + C\dot{\gamma} + 2Fa_2D_r], \quad (4)$$

where P is the total hydrostatic pressure and $D_r$ is the rotary diffusivity due to Brownian motion. For slender particles (i.e., $r \gg 1$), the fibers are large enough that the Brownian motion can be omitted [Chaouche and Koch (2001)]. For semiconcentrated fiber suspensions, Dinh and Armstrong (1984) have calculated the forces induced by the fluid on cylinders of an infinite aspect ratio and derived the particle motion and rotation. Using the cell model approach developed by Batchelor (1970), they expressed the stress contribution of the fibers by averaging the fiber contribution with the orientation distribution function. Equation (4) reduces to that of the Dinh and Armstrong (1984) result if A is the only nonzero coefficient.

More recently, fiber interactions have been modeled by direct simulation methods. These techniques have been used for rigid fibers by Yamamoto and Matsuoka (1996) and Ausias et al. (2006) who took also into account a fiber length distribution. Schmid et al. (2000) and Joung et al. (2001) did similar investigations for flexible fibers. Sundararajakumar and Koch (1997) investigated the effects of direct fiber contacts by neglecting the hydrodynamic interactions. It was observed an increase of fiber orientation resulting in an enhancement of the shear viscosity and the normal stress differences. They also concluded that the hydrodynamic interactions remain the basic aspect of the structure in dilute and semidilute fiber suspensions. Although of considerable interest, these direct numerical simulations cannot be easily introduced into commercial software packages owing to their large computational time.

Finally, Sandstrom (1993) and Sandstrom and Tucker (1993) have presented a theory for concentrated fiber suspensions with strong fiber interactions. The interfiber contact forces are modeled as linear lubrication forces between particles and are weighted by the probability distribution function. Nonlinear formulations have been proposed by Servais et al. (1999a, 1999b) and Djalili-Moghaddam and Toll (2005). Le Corre et al. (2005) used a numerical analysis based on the homogenization of discrete and periodic structures to investigate the relationships between the microstructure and the macroscopic rheological properties of concentrated suspensions of fiber bundles having an average planar orientation. More developments have been achieved by Djalili-Moghaddam and Toll (2005) who considered the total stress as the sum of the contributions of the matrix, the hydrodynamic forces acting on the particles and the interaction forces between fibers.

In this paper we present a rheological model for semiconcentrated fiber suspensions, in which interactions between fibers are taken into account. In the case of dilute suspensions the interactions will become negligible and the model will reduce to that of Dinh and Armstrong (1984). The approach uses the fiber orientation distribution function from which new conformation tensors, called interaction tensors, are constructed. An evolution
<figure><img src="image_1.png" /><figcaption>FIG. 1. Schematic representation of two interacting fibers. Superscript α represents the test fiber, superscript β refers to the neighboring fiber.</figcaption></figure>

equation for the orientation is derived, and a constitutive law is developed by adding a fiber-fiber contribution to the total stress equation proposed by Dinh and Armstrong (1984). Closure approximations required to evaluate the fourth-order interaction tensor are also proposed. The new model contains a few adjustable parameters that are determined by fitting experimental data for short glass fiber-filled reinforced polybutene.

## II. MODEL FORMULATION

### A. Definitions

The test fiber is denoted by the superscript $\alpha$, whereas the superscript $\beta$ refers to the neighboring fiber (Fig. 1). The location of the mass center of the test fiber is $\mathbf{r}_C^{\alpha}$, and $s^{\alpha}$ is an arc length measured along the fiber axis with $s^{\alpha}=0$ at the center. The unit vector $\mathbf{p}^{\alpha}$ along the major axis of the fiber describes its orientation. The vector $\mathbf{r}_C^{\alpha}$ is expressed with respect to a fixed coordinate system $\mathbf{R}$ (the velocity of the coordinate origin is set to zero) and a particle coordinate system $\mathbf{R}^{\alpha}$ (where the 1-direction is taken collinear to $\mathbf{p}^{\alpha}$) on the fiber element. Similar definitions are used for the neighboring fiber using the superscript $\beta$. The surrounding fluid angular velocity is defined by $\Omega^{\infty} = \frac{1}{2} \nabla \times \mathbf{v}^{\infty}$.

### B. Hypotheses

(1) All fibers are considered as straight and rigid cylinders of uniform length $L$, constant diameter $D$, and large aspect ratio $r \gg 1$ (absence of Brownian motion).

(2) The bulk flow is assumed to be homogeneous, whereas the velocity gradient is constant over the fiber length ($\nabla \mathbf{v}^{\infty} = \mathbf{\kappa}^{\dagger}$), so the fluid velocity is given by $\mathbf{v}^{\infty} = \mathbf{\kappa} \cdot \mathbf{r}_C^{\alpha}$.

(3) The suspended fibers are assumed to be distributed uniformly, i.e., there is no concentration gradient. Therefore, for homogeneous flows the probability distribution functions $\psi_{p^{\alpha}}$ and $\psi_{p^{\beta}}$ are independent on the spatial location of the fibers and $n$ becomes a constant.

(4) The translational and rotational inertia of particles are neglected. Furthermore, gravitational effects are negligible.
(5) The suspending fluid is assumed to be incompressible and a continuous Newtonian medium. The fluid motion is determined by the Stokes equation (negligible translational and rotational fluid inertia).

(6) The particle-particle interaction force is linear with respect to the relative velocity at the contact point (hydrodynamic friction).

(7) Affine motion and rotation are then assumed.

### C. Equations for fiber motion

From Kim and Karrila (1991), the time evolution for the change of orientation of a test prolate spheroid $\alpha$ yields

$$ \dot{\mathbf{p}}_{H}^{\alpha}=\mathbf{\omega}^{\alpha} \times \mathbf{p}^{\alpha}=(\mathbf{\Omega}^{\infty}+\mathbf{C}^{-1} \cdot \mathbf{H}: \dot{\mathbf{\gamma}}) \times \mathbf{p}^{\alpha}=-\frac{1}{2} \mathbf{\omega} \cdot \mathbf{p}^{\alpha}+\frac{1}{2} \frac{Y_{H}}{Y_{C}}\left(\dot{\mathbf{\gamma}} \cdot \mathbf{p}^{\alpha}-\dot{\mathbf{\gamma}}: \mathbf{p}^{\alpha} \mathbf{p}^{\alpha}\right), \quad (5) $$

where $\mathbf{\omega}^{\alpha}$ is the angular velocity of fiber $\alpha$. Expressed in $\mathbf{R}$, the resistance tensors are, $\mathbf{C}=\eta_{0} L^{3}[Y_{C} \mathbf{\delta}+(X_{C}-Y_{C}) \mathbf{p}^{\alpha} \mathbf{p}^{\alpha}]$ and $\mathbf{H}=-\frac{1}{2} \eta_{0} L^{3} Y_{H}(\mathbf{\epsilon} \cdot \mathbf{p}^{\alpha}) \mathbf{p}^{\alpha}$, respectively, where $\mathbf{\epsilon}$ is the third order permutation tensor and $X_C, Y_C$, and $Y_H$ are resistance functions. For the dilute regime, we find that the rotational drag coefficients are $X_C=4\pi/r^2$ and $Y_C=1/6$. Expression (5) is equivalent to the Jeffery (1922) equation, when $Y_H/Y_C$ is the form factor [Kim and Karrila (1991)]. In the limit of infinite aspect ratio (i.e., the form factor, $\lambda=1$), the slender body theory is recovered, which reduces the Jeffery equation to

$$ \dot{\mathbf{p}}_{H}^{\alpha}=-\frac{1}{2} \mathbf{\omega} \cdot \mathbf{p}^{\alpha}+\frac{1}{2}\left(\dot{\mathbf{\gamma}} \cdot \mathbf{p}^{\alpha}-\dot{\mathbf{\gamma}}: \mathbf{p}^{\alpha} \mathbf{p}^{\alpha}\right)=\mathbf{\kappa} \cdot \mathbf{p}^{\alpha}-\mathbf{\kappa}: \mathbf{p}^{\alpha} \mathbf{p}^{\alpha}. \quad (6) $$

For the case of fiber motion with interactions, we assume that the test fiber $\alpha$ suspended in a Newtonian fluid is subjected to an external force caused by a neighboring interacting fiber $\beta$. Using the previous development, the equation for the change of orientation for any fiber $\alpha$ becomes

$$ \dot{\mathbf{p}}^{\alpha}=\mathbf{\omega}^{\alpha} \times \mathbf{p}^{\alpha}=[\mathbf{\Omega}^{\infty}+\mathbf{C}^{-1} \cdot\left(\mathbf{H}: \dot{\mathbf{\gamma}}+\mathbf{T}_{1}\right)] \times \mathbf{p}^{\alpha}, \quad (7) $$

where $\mathbf{T}_1$ represents the global torque produced by all neighboring fibers, which act on the test fiber $\alpha$. The three terms on the right-hand side of Eq. (7) refer to the rotations due to the rigid body, the fluid deformation, and the interaction torque, respectively. Although interactions between two particles are very complex, some specific models exist when dealing with suspensions of prolate spheroids in Stokes flows [Claeys and Brady (1993)]. In this work, a less tedious and a more phenomenological for cylindrical shape model has been used. To determine $\mathbf{T}_1$, the differential interaction moment $d\mathbf{T}_1$ has to be used, which is given by $d\mathbf{T}_1=s^{\alpha}\mathbf{p}^{\alpha} \times d\mathbf{f}_1$. It is simply the lever arm times the applied differential force $d\mathbf{f}_1$. The general form for the force due to short range fiber-fiber interactions has been proposed by several authors [Sandstrom (1993); Toll and Manson (1994); Gibson and Toll (1999); Djalili-Moghaddam and Toll (2005)]. If a linear lubrication is assumed, it gives

$$ \mathbf{f}_1=D \eta_0 k \Delta \dot{\mathbf{r}}, \quad (8) $$

where $k$ is a dimensionless geometric factor and $\Delta\dot{\mathbf{r}}$ is the time derivative of the shortest distance between two fibers $\Delta\mathbf{r}$. Contrary to Djalili-Moghaddam and Toll (2005), who deal with short-range interactions (i.e., $|\Delta\mathbf{r}| \ll s^{\alpha}$), no term is neglected in this work to express the relative velocity $\Delta\dot{\mathbf{r}}$ at the interaction point. Then, contacts close to fiber mass center are taken into account. Therefore, $\Delta\dot{\mathbf{r}}$ is written in a fixed coordinate system as
$$ 
\begin{aligned}
\Delta \dot{\mathbf{r}} ={}& \dot{\mathbf{r}}_{C}^{\beta} + s^{\beta} \dot{\mathbf{p}}^{\beta} - \dot{\mathbf{r}}_{C}^{\alpha} - s^{\alpha} \dot{\mathbf{p}}^{\alpha} = \dot{\mathbf{r}}_{C}^{\beta} - \dot{\mathbf{r}}_{C}^{\alpha} + \frac{s^{\beta}}{2} (-\boldsymbol{\omega} \cdot \mathbf{p}^{\beta} + \dot{\boldsymbol{\gamma}} \cdot \mathbf{p}^{\beta} - \dot{\boldsymbol{\gamma}} : \mathbf{p}^{\beta} \mathbf{p}^{\beta} \mathbf{p}^{\beta}) \\
& - \frac{s^{\alpha}}{2} (-\boldsymbol{\omega} \cdot \mathbf{p}^{\alpha} + \dot{\boldsymbol{\gamma}} \cdot \mathbf{p}^{\alpha} - \dot{\boldsymbol{\gamma}} : \mathbf{p}^{\alpha} \mathbf{p}^{\alpha}).
\end{aligned}
\quad (9) $$

To determine the fiber rotation only, $\Delta \dot{\mathbf{r}}$ has to be expressed in a corotational frame (this frame moves and rotates with a fluid element) in order to satisfy the requirement of material objectivity. Rigid-body rotation is already taken into account by the term $\Omega^{\infty} \times \mathbf{p}^{\alpha}$ in Eq. (7). Therefore, in this corotational frame (denoted by the superscript crf), $\Delta \dot{\mathbf{r}}$ reduces to the following form:

$$ \Delta \dot{\mathbf{r}}^{\text{crf}} = \frac{s^{\beta}}{2} (\dot{\boldsymbol{\gamma}} \cdot \mathbf{p}^{\beta} - \dot{\boldsymbol{\gamma}} : \mathbf{p}^{\beta} \mathbf{p}^{\beta}) - \frac{s^{\alpha}}{2} (\dot{\boldsymbol{\gamma}} \cdot \mathbf{p}^{\alpha} - \dot{\boldsymbol{\gamma}} : \mathbf{p}^{\alpha} \mathbf{p}^{\alpha}). \quad (10) $$

$\Delta \dot{\mathbf{r}}^{\text{crf}}$ becomes a simple function of the rate-of-deformation tensor. It means that only the fluid deformation modifies the distance $\Delta \mathbf{r}$ between two fibers at the interaction point. The normal component of the relative velocity $\Delta \dot{\mathbf{r}}^{\text{crf}}$ with respect to the surface formed by two contacting fibers is $\Delta \dot{\mathbf{r}}^{\text{crf}} \cdot (\mathbf{p}^{\beta} \times \mathbf{p}^{\alpha})$, and this term is not null compared to the relative velocity relation used by Djalili-Moghaddam and Toll (2005).

The expression for $\Delta \dot{\mathbf{r}}^{\text{crf}}$ is then inserted into Eq. (8) (interaction force) to determine the fiber rotation. In the next step we consider that the short-range interaction force is function of the fiber contact probability. As the average interfiber spacing increases, one would expect the effect of interaction to decrease and vice versa. The interfiber spacing is a physical parameter that directly indicates the available space for fiber rotation. The interacting fiber $\beta$ modifies the test fiber $\alpha$ rotation through the force weighted by the probability that the neighboring fiber interacts with the test fiber. This quantity is proportional to the probability $\psi_{\mathbf{p}^{\beta}}$ to find a fiber with an orientation $\mathbf{p}^{\beta}$. The orientation of the neighboring fiber is independent of that of the test fiber, which means that the orientations of neighboring fibers are not correlated. We used the following classical assumptions: a lubricated contact occurs if the distance between the nearest points of the two fibers is shorter than the fiber diameter $D$ and there is no excluded volume effect. Hence, the model is restricted to semiconcentrated suspensions. The positions along the test fiber and the interacting fiber of this nearest point are $s^{\alpha}$ and $s^{\beta}$, respectively (see Fig. 1). Doi and Edwards (1978), and then Sandstrom (1993) have proposed an evaluation of this probability as

$$ P_C = 2nD|\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}|\psi_{\mathbf{p}^{\beta}}d\mathbf{p}^{\beta}ds^{\alpha}ds^{\beta}. \quad (11) $$

The interaction force [Eq. (8)] is weighted by this contact probability $P_C$ to give the elementary interaction force

$$ d\mathbf{f}_1 = 2k\eta_0 n D^2 \Delta \dot{\mathbf{r}}^{\text{crf}} |\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}| \psi_{\mathbf{p}^{\beta}} d\mathbf{p}^{\beta} ds^{\alpha} ds^{\beta}. \quad (12) $$

Then, the above expression is used to express the global torque $\mathbf{T}_1$. Integration over the neighboring and test fiber lengths gives all the possible interaction positions. It implies that all the contact point positions (described by the arc lengths $s^{\alpha}$ and $s^{\beta}$) will have the same probability to occur. An extra integration over all directions of $\mathbf{p}^{\beta}$ takes into account all the possible orientations of the neighboring fiber. With these considerations, the global torque $\mathbf{T}_1$ is evaluated to result in the expression of the orientation time evolution for the test fiber $\alpha$,
$$ \dot{\mathbf{p}}^{\alpha} = \dot{\mathbf{p}}_{H}^{\alpha} - \phi \tilde{M} \int_{\mathbf{p}^{\beta}} |\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}| \left( \frac{1}{2} \dot{\gamma} \cdot \mathbf{p}^{\alpha} - \frac{1}{2} \dot{\gamma} : \mathbf{p}^{\alpha} \mathbf{p}^{\alpha} \right) \psi_{\mathbf{p}^{\beta}} d\mathbf{p}^{\beta}, \quad (13) $$

where the dimensionless scalar $\tilde{M}=2k/3\pi Y_c$. The term in brackets represents the fiber rotation induced by fluid deformation only. Since interactions modify also fiber rotation, a perturbation velocity is introduced similarly as done by Folgar and Tucker (1984), which is proportional to the effective deformation rate $|\dot{\gamma}|$. This diffusional term is not due to the Brownian motion, but represents the rotational fluctuations undergoing by the test fiber, when the suspension is flowing. Therefore, the test fiber evolution can be expressed as

$$ \dot{\mathbf{p}}^{\alpha} = \dot{\mathbf{p}}_{H}^{\alpha} - \phi \tilde{M} \int_{\mathbf{p}^{\beta}} |\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}| \left( \frac{1}{2} \dot{\gamma} \cdot \mathbf{p}^{\alpha} - \frac{1}{2} \dot{\gamma} : \mathbf{p}^{\alpha} \mathbf{p}^{\alpha} + q |\dot{\gamma}| \frac{\partial \ln \psi_{\mathbf{p}^{\alpha}}}{\partial \mathbf{p}^{\alpha}} \right) \psi_{\mathbf{p}^{\beta}} d\mathbf{p}^{\beta}, \quad (14) $$

where $q$ is a dimensionless interaction coefficient. We notice that when the fiber volume fraction $\phi$ is small, Eq. (14) reduces to the time evolution equation for slender particles [Eq. (6)].

### D. Time evolution of the second order orientation tensor

From Bird et al. (1987), the probability distribution function may be regarded as a convected quantity. Then, the dynamic change of the second orientation tensor is derived from Eq. (14) as

$$ \frac{D\mathbf{a}_2}{Dt} = -\frac{1}{2}(\boldsymbol{\omega} \cdot \mathbf{a}_2 - \mathbf{a}_2 \cdot \boldsymbol{\omega}) + \frac{1}{2}(\dot{\boldsymbol{\gamma}} \cdot \mathbf{a}_2 + \mathbf{a}_2 \cdot \dot{\boldsymbol{\gamma}} - 2\dot{\boldsymbol{\gamma}}:\mathbf{a}_4) - \frac{1}{2}\phi \tilde{M}(\dot{\boldsymbol{\gamma}} \cdot \mathbf{b}_2 + \mathbf{b}_2 \cdot \dot{\boldsymbol{\gamma}} - 2\dot{\boldsymbol{\gamma}}:\mathbf{b}_4) \\ - \phi \tilde{M} q |\dot{\boldsymbol{\gamma}}| \int_{\mathbf{p}^{\alpha}} \frac{\partial \psi_{\mathbf{p}^{\alpha}}}{\partial \mathbf{p}^{\alpha}} \int_{\mathbf{p}^{\beta}} |\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}| \psi_{\mathbf{p}^{\beta}} d\mathbf{p}^{\beta} \cdot \frac{\partial \mathbf{p}^{\alpha}}{\partial \mathbf{p}^{\alpha}} d\mathbf{p}^{\alpha}. \quad (15) $$

Some already known tensors appear, such as the second- and fourth-order orientation tensors $\mathbf{a}_2$ and $\mathbf{a}_4$ [Advani and Tucker (1987)]. New tensors are also introduced, the second- and fourth-order tensors $\mathbf{b}_2$ and $\mathbf{b}_4$, called interaction tensors. They are defined by forming the dyadic products of the vector $\mathbf{p}^{\alpha}$, weighted by the potential $|\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}|$, and then twice integrating the product with respect to the distribution function over all possible directions [Djalili-Moghaddam and Toll (2005)]. Because the distribution function is even, the odd-order integrals are zero, so only the even-order tensors are of interest

$$ \mathbf{b}_2 = \int_{\mathbf{p}^{\alpha}} \int_{\mathbf{p}^{\beta}} \mathbf{p}^{\alpha} \mathbf{p}^{\alpha} |\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}| \psi_{p\beta} \psi_{p\alpha} d\mathbf{p}^{\beta} d\mathbf{p}^{\alpha}, \quad (16) $$

$$ \mathbf{b}_4 = \int_{\mathbf{p}^{\alpha}} \int_{\mathbf{p}^{\beta}} \mathbf{p}^{\alpha} \mathbf{p}^{\alpha} \mathbf{p}^{\alpha} \mathbf{p}^{\beta} |\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}| \psi_{p\beta} \psi_{p\alpha} d\mathbf{p}^{\beta} d\mathbf{p}^{\alpha}. \quad (17) $$

The interaction tensors are completely symmetric. Moreover, all the terms in Eq. (15) have been averaged easily except the last term on the right-hand side, which requires an approximate solution. To achieve a simplest form, we follow the approach done for the FENE-P model by Peterlin (1955), Bird et al. (1987), and Carreau et al. (1997) who consider the average macromolecule length as a characteristic length. In our case, the preaveraging of the scalar quantity $\int_{\mathbf{p}^{\beta}} |\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}| \psi_{\mathbf{p}^{\beta}} d\mathbf{p}^{\beta}$ with respect to the distribution
<table><caption>TABLE I. Calculated values of *f* and *g* for different fiber orientations.</caption><thead><tr><th>Orientation</th><th>f</th><th>g</th></tr></thead><tbody><tr><td>Unidirectional</td><td>0</td><td>1</td></tr><tr><td>2D random</td><td>2/π</td><td>2π</td></tr><tr><td>2D (1/2 → 1; 1/2 → 2)<sup>a</sup></td><td>1/2</td><td>1/2</td></tr><tr><td>3D random</td><td>π/4</td><td>1/2</td></tr><tr><td>3D (1/3 → 1; 1/3 → 2; 1/3 → 3)<sup>b</sup></td><td>2/3</td><td>1/3</td></tr></tbody></table>

<sup>a</sup>1/2 → 1 ≡ 50% of fibers aligned in the 1-direction.

<sup>b</sup>1/3 → 1 ≡ 1/3 of fibers aligned in the 1-direction.

$\psi_{\mathbf{p}^\alpha}$ indicates that the average contact number is representative of all the contact points. Hence, we obtain the following scalar function:

$$f = \text{trace } \mathbf{b}_2 = \int_{\mathbf{p}^\alpha} \int_{\mathbf{p}^\beta} |\mathbf{p}^\alpha \times \mathbf{p}^\beta| \psi_{\mathbf{p}^\beta} \psi_{\mathbf{p}^\alpha} d\mathbf{p}^\beta d\mathbf{p}^\alpha. \quad (18)$$

Toll (1998) referred to *f* as the average number of contacts per fiber. Moreover, he defined a second scalar invariant *g* of the fiber orientation distribution, defined as

$$g = \int_{\mathbf{p}^\alpha} \int_{\mathbf{p}^\beta} |\mathbf{p}^\alpha \cdot \mathbf{p}^\beta| \psi_{\mathbf{p}^\beta} \psi_{\mathbf{p}^\alpha} d\mathbf{p}^\beta d\mathbf{p}^\alpha. \quad (19)$$

Some generic values of *f* and *g* for different fiber orientations are listed in Table I. According to Toll (1998), the scalar quantities *f* and *g* are related to the average number of lubricated contacts per fiber $N_C$ as

$$N_C = \frac{8}{\pi} r \phi f + 4 \phi (g + 1). \quad (20)$$

For slender fibers, the average number of contacts per fiber reduces to

$$N_C = \frac{8}{\pi} r \phi f. \quad (21)$$

Using the preaveraging, the equation of change for $\mathbf{a}_2$ becomes

$$\frac{D\mathbf{a}_2}{Dt} = -\frac{1}{2}(\boldsymbol{\omega} \cdot \mathbf{a}_2 - \mathbf{a}_2 \cdot \boldsymbol{\omega}) + \frac{1}{2}(\dot{\boldsymbol{\gamma}} \cdot \mathbf{a}_2 + \mathbf{a}_2 \cdot \dot{\boldsymbol{\gamma}} - 2\dot{\boldsymbol{\gamma}}:\mathbf{a}_4) - \frac{1}{2}\phi \tilde{M}(\dot{\boldsymbol{\gamma}} \cdot \mathbf{b}_2 + \mathbf{b}_2 \cdot \dot{\boldsymbol{\gamma}} - 2\dot{\boldsymbol{\gamma}}:\mathbf{b}_4) + 2f\phi \tilde{M}q|\dot{\boldsymbol{\gamma}}|(\boldsymbol{\delta} - 3\mathbf{a}_2), \quad (22)$$

where the last term is a diffusion term. If we refer to Eq. (14), which gives the average angular velocity of a particle in a concentrated suspension, the variance is equal to $2f\phi \tilde{M}q|\dot{\boldsymbol{\gamma}}|$ compared to $2C_1|\dot{\boldsymbol{\gamma}}|$ for the Folgar and Tucker (1984) model. This quantity evolves with time: it is maximum for a random orientation state and tends to zero for perfectly aligned fibers (see *f* values in Table I). The variance in this model is a function of the number of fiber contacts, not accounted for in the Folgar and Tucker model. As demonstrated in Appendix A, Eq. (22) has the same structure than the reduced-strain closure model, which was recently developed to slow down the fiber orientation kinetics [Tucker *et al.* (2007) and Phelps *et al.* (2008)].

In Appendix B, it is proven that the average mass center of a fiber moves affinely with the bulk flow.
### III. STRESS DETERMINATION

Based on the work of Dinh and Armstrong (1984) for slender fiber suspensions, we express the total stress tensor as

$$ \boldsymbol{\sigma} = -P\boldsymbol{\delta} + \eta_0\dot{\boldsymbol{\gamma}} + \eta_0\phi\frac{r^2}{6\pi}X_A\dot{\boldsymbol{\gamma}}:\boldsymbol{a}_4, \quad (23) $$

where $X_A$ is the parallel drag coefficient of the fibers. It depends on the nature of the fiber-matrix contact and is taken as an adjustable parameter in this work. For semiconcentrated suspensions, there is an additional stress contribution due to the contact forces acting on the fibers. From the previous development, the deviatoric stress $\boldsymbol{\tau}$ becomes

$$ \boldsymbol{\tau} = \eta_0 \dot{\boldsymbol{\gamma}} + \eta_0 \phi \frac{r^2}{6\pi} X_A \dot{\boldsymbol{\gamma}}:\boldsymbol{a}_4 + n \iint_{\boldsymbol{p}^\alpha} s^\alpha \boldsymbol{p}^\alpha d\boldsymbol{f}_1 \psi_{\boldsymbol{p}^\alpha} d\boldsymbol{p}^\alpha, \quad (24) $$

where the stress contribution comes from the Newtonian matrix, the hydrodynamic forces (same as previously), and the fiber interactions [Toll and Manson (1994)], respectively [hypothesis (7)]. Combining Eqs. (9), (12), and (24) yields the stress induced by fiber interactions

$$ \boldsymbol{\tau}_1 = \eta_0 \phi^2 \frac{4r^2}{3\pi^2} k \dot{\boldsymbol{\gamma}}:\boldsymbol{b}_4. \quad (25) $$

We note that the fluid velocity at any contact point is close to that of the fluid velocity at a neighboring contact point. Hence, using the assumption that particles do not perturb the fluid flow, the velocity at the contact point is $\dot{\boldsymbol{r}}_C^\beta + s^\beta \boldsymbol{\kappa} \cdot \boldsymbol{p}^\beta \approx \dot{\boldsymbol{r}}_C^\alpha + s^\alpha \boldsymbol{\kappa} \cdot \boldsymbol{p}^\alpha$. Finally, the total stress tensor including fiber-fiber interactions is given by

$$ \boldsymbol{\sigma} = -P\boldsymbol{\delta} + \eta_0 \dot{\boldsymbol{\gamma}} + \eta_0 \phi \frac{r^2}{6\pi} X_A \dot{\boldsymbol{\gamma}}:\boldsymbol{a}_4 + \eta_0 \phi^2 \frac{4r^2}{3\pi^2} k \dot{\boldsymbol{\gamma}}:\boldsymbol{b}_4. \quad (26) $$

Expression (26) reduces to that proposed by Djalili-Moghaddam and Toll (2005) for short-range interactions, if the slender-body coefficients (used here) are replaced by the Shaqfeh-Fredrickson coefficients [Shaqfeh and Fredrickson (1990)]. Furthermore, in the work of Djalili-Moghaddam and Toll (2005) the interaction distance is arbitrary, whereas in this work it is taken as the fiber diameter, as already done by Sandstrom (1993). Again, we note that for very low fiber volume fraction, the term in $\phi^2$ becomes insignificant and the total stress [Eq. (26)] reduces to the Dinh and Armstrong (1984) result [Eq. (23)].

### IV. APPROXIMATIONS

Equations (22) and (26) do not represent a complete set of rheological equations yet. To achieve this, we need to state some approximations. The first one has to relate the fourth-order orientation tensor to the second order one. Many closure approximations have been proposed in the literature, but will not be reported further in this paper [see Chung and Kwon (2001) and Sepehr et al. (2004b) for more details]. The ORW3 closure approximation (orthotropic fitted closure approximation for wide $C_I$ values with third-order polynomial expansions) developed by Chung and Kwon (2001) is used to evaluate the fourth-order orientation tensor. This closure has been chosen as the diffusion coefficient in Eq. (22) evolves with time and it has been shown to give reliable results for a wide range of $C_I$ [Chung and Kwon (2001, 2002)]. In the next development, a relation-
ship between the second-order orientation tensor and the second-order interaction tensor is proposed. Finally, we suggest closure approximations for the fourth-order interaction tensor.

## A. Onsager and Maier-Saupe potentials

First, we note that the scalar quantity $|\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}|$ known as the Onsager potential [Onsager (1949)] cannot be expressed in a simple form without the use of a closure approximation. Nevertheless, it is possible to relate the interaction tensor $\mathbf{b}_2$ to the orientation tensors (second and fourth) by correctly approximating the Onsager potential. It is often replaced by the Maier-Saupe potential [Maier and Saupe (1959)] expressed as $|\mathbf{p}^{\alpha} \times \mathbf{p}^{\beta}|^2$. These two potentials are obviously different, but they share the essential properties that their maximum are reached when the fibers are orthogonal and then decrease when the fibers orient [Grmela and Dlugogorski (1996)]. By performing some straightforward algebraic calculations (see Appendix C), we obtain

$$ \mathbf{b}_2 = \frac{3\pi}{8}(\mathbf{a}_2 - \mathbf{a}_4:\mathbf{a}_2), \quad (27) $$

where the numerical prefactor $3\pi/8$ is introduced to reach the exact value of $f=\pi/4$ for an isotropic fiber orientation, as already done by Dhont and Briels (2003) who evaluated the Onsager potential using a Ginzburg-Landau expansion.

## B. Interaction closure approximations

As initial simple approximations we propose two closures, one based on the quadratic formulation and the second one formulated with a linear combination. The quadratic closure, called IQUA, has the following form:

$$ b_{ijkl}^{Q} = \frac{1}{f} b_{ij} b_{kl}, \quad (28) $$

whereas the linear closure (ILIN) can be written as

$$ b_{ijkl}^{L} = -\frac{f}{35}(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \frac{1}{7}(b_{ij}\delta_{kl} + b_{ik}\delta_{jl} + b_{il}\delta_{jk} + b_{kl}\delta_{ij} + b_{jl}\delta_{ik} + b_{jk}\delta_{il}). \quad (29) $$

It was shown by Sepehr et al. (2004a) that linear orientation closures produced oscillations, while quadratic orientation approximations were not efficient to predict the correct behavior for reverse flows. Results obtained with the quadratic form (IQUA) seem to be more relevant than the ones obtained with the linear interaction closure (ILIN). Hence, for simplicity, we chose the IQUA approximation for the interaction tensor in the following comparison with experimental results. Other closure approximations will be examined in another investigation.

# V. PREDICTIONS AND COMPARISON WITH EXPERIMENTAL DATA

## A. Model predictions

The present model, formed by Eqs. (22) and (26)-(28), has been tested in simple shear flow for forward and reverse stress growth experiments. The deformation is calculated as $\gamma=\dot{\gamma}t$, where the applied shear rate is set at $1 \text{ s}^{-1}$. All predictions presented in this section were obtained for the following conditions: $\eta_0=1 \text{ Pa} \cdot \text{s}$, $\phi=10\%$ and $r=20$. The closure
<figure><img src="image_2.png" /><figcaption>FIG. 2. Predicted stress growth functions for forward flow (CW); ORW3 and IQUA closure approximations, $\dot{\gamma}=1$ s⁻¹ and $q=0.5$: (a) transient viscosity; (b) transient primary normal stress differences as functions of strain.</figcaption></figure>

approximation for the fourth-order interaction tensors is the IQUA [Eq. (28)]. The only adjustable parameters are the dimensionless numbers $k$ and $q$, with $k$ taking into account the intensity of the friction force and $q$ proportional to the diffusion term. The resistance functions $X_A$ and $Y_C$ are set to 1 and 1/6, respectively, from the dilute regime. Results from the Folgar–Tucker–Dinh–Armstrong (FTDA) model are also compared. The FTDA model consists of Eqs. (1) and (23), where the interaction coefficient $C_I$ is set to 0.01. For both models, the initial fiber orientation is assumed to be three-dimensional (3D) random (isotropic) and the closure approximation for the fourth-order orientation tensor is the ORW3 [Chung and Kwon (2001)]. The system of equations of each model is solved using the MATLAB® software.

Figures 2 and 3 report the predictions for the shear stress growth coefficient $\eta^+$ and the first normal stress growth function $N_1^+$ in forward [clockwise (CW) direction] and reverse [counterclock-wise (CCW) direction] flows, respectively, for different values of $k$. In comparison to the FTDA predictions, the results show that by increasing $k$ the steady-state values are reached more rapidly. In addition, the relative amplitude of the over-
<figure><img src="image_3.png" /><figcaption>FIG. 3. Predicted stress growth functions for reverse flow (CCW); ORW3 and IQUA closure approximations, $|\dot{\gamma}|=1$ s$^{-1}$ and $q=0.5$: (a) transient viscosity; (b) transient primary normal stress differences as functions of strain.</figcaption></figure>

shoots for the stress growth functions in the forward direction [Figs. 2(a) and 2(b)] depends on the $k$ value, especially for the first normal stress growth function. We note some oscillations in the predicted functions. For the reverse direction [Figs. 3(a) and 3(b)], the overshoots are predicted at lower strain values as $k$ increases. It is a reasonable result since fewer fibers are aligned in the flow direction for large $k$. All of these observations confirm that fiber interactions prevent fiber orientation, which is physically realistic. As in Figs. 2(a) and 2(b), we note oscillations in the predicted functions for the reverse (ccw) flow: the local undershoots followed by a slight overshoot, observed at strains larger than 10 are probably due to the quadratic closure approximation used for the fourth-order interaction tensor.

Effects of different values of $q$ on $\eta^+$ and $N_1^+$ are presented in Figs. 4 and 5 for forward and reverse flow experiments, respectively. In the forward direction [Figs. 4(a) and 4(b)], $q$ does not affect the strain at which the overshoots of the stress growth functions are predicted. In addition, as $q$ increases the maximum value of the stress growth coefficient $\eta^+$ is unchanged, whereas the maximum of the first normal stress growth function $N_1^+$
<figure><img src="image_4.png" /><figcaption>FIG. 4. Predicted stress growth functions for forward flow (CW); ORW3 and IQUA closure approximations, $\dot{\gamma}=1$ s$^{-1}$ and $k=0.3$: (a) transient viscosity; (b) transient primary normal stress differences as functions of strain.</figcaption></figure>

slightly decreases. However, the relative amplitude of the overshoots with respect to the steady-state values is a strong function of $q$ and the steady values become larger as $q$ is increased. In analogy with results of Figs. 3(a) and 3(b), the overshoots for reverse flows [Figs. 5(a) and 5(b)] occur at lower strain values as $q$ increases. Similar results have been reported by Sepehr et al. (2004b) using the Folgar–Tucker–Lipcomb model and different values of the interaction coefficient $C_I$. Since parameter $q$ is related to diffusion, it is reasonable to expect that $q$ and $C_I$ produce similar effects.

Key components of $\mathbf{a}_2$ and the trace of $\mathbf{b}_2$, $f$, have been also calculated for forward (Fig. 6) and reverse (Fig. 7) stress growth experiments, both for $k=0.3$ and $q=0.5$. Figure 6 indicates that for an initial 3D random state most of the fibers orient in the shear direction. The scalar $f$ depicts a decrease from $f=\pi/4$ to a steady value of $f=0.4588$, with a small undershoot at a strain close to 10. For the reverse flow (Fig. 7), fiber tumbling is predicted by the model, as illustrated by the $\mathbf{a}_2$ components. Also, $f$ exhibits an overshoot before reaching the same steady state of 0.4588 indicating that fiber contacts increase during tumbling.
<figure><img src="image_5.png" /><figcaption>FIG. 5. Predicted stress growth functions for reverse flow (CCW); ORW3 and IQUA closure approximations, $|\dot{\gamma}|=1$ s⁻¹ and $k=0.3$: (a) transient viscosity; (b) transient primary normal stress differences as functions of strain.</figcaption></figure>

## B. Comparison with experimental data

Experiments have been performed on suspensions of short glass fibers in a Newtonian polybutene matrix (PB Indopol H100, η₀=35.1 Pa·s at 20 °C). The fibers were obtained from the calcination at 600 °C for 40 min of a 30 wt % filled polypropylene provided by Basell. The average fiber length *L* was close to 320 µm and its diameter *D*=16 µm. Four different weight fractions of fibers were used [10% (PB10), 15% (PB15), 20% (PB20), and 25% (PB25), corresponding to volume fractions, *φ*, of 3.8%, 5.8%, 8.1% and 10.5%, respectively]. Sample PB10 was in the semidilute regime whereas PB15, PB20, and PB25 were in the concentrated one. The composites were prepared manually by gently mixing the fibers in the matrix in order to avoid fiber breakage. Measurements were performed with an Anton Paar Physica MCR501 rheometer using a parallel plate geometry with a 25 mm radius. The gap was fixed at 1.5 mm, which allowed fibers to rotate freely at least in the shear plane. Stress growth experiments were carried out at a shear rate of 10 s⁻¹. The torque and the normal force exerted by the material on the disks are, respectively, given by [Carreau *et al.* (1997)];
<figure><img src="image_6.png" /><figcaption>FIG. 6. Predicted variations for key components of $\mathbf{a}_2$ ($a_{11}$, $a_{22}$, and $a_{12}$) and for the scalar $f$ as functions of deformation for forward flow (ORW3 and IQUA closure approximations, $\dot{\gamma}=1$ s$^{-1}$, $X_A=1$, $k=0.3$, and $q=0.5$).</figcaption></figure>

$$T = \int_{0}^{2\pi} \int_{0}^{R} \sigma_{\theta z}(r) r^2 dr d\theta, \qquad (30)$$

$$F = \int_{0}^{2\pi} \int_{0}^{R} \sigma_{zz}(r) r dr d\theta. \qquad (31)$$

These quantities are dependent on the strain $\gamma$ and different values are obtained at different radius $r$. To compare the experimental data with the model predictions, the predicted rheological functions have to be evaluated under the same conditions, i.e., same deformation. Using straightforward calculations, the model predictions are calculated as,

<figure><img src="image_7.png" /><figcaption>FIG. 7. Predicted variations for key components of $\mathbf{a}_2$ ($a_{11}$, $a_{22}$, and $a_{12}$) and for the scalar $f$ as functions of deformation for reverse flow (ORW3 and IQUA closure approximations, $|\dot{\gamma}|=1$ s$^{-1}$, $X_A=1$, $k=0.3$, and $q=0.5$).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 8. Predicted reduced steady shear viscosity (a) and normal stress differences (b) compared to experimental data as functions of fiber volume fraction (ORW3 and IQUA closure approximations). Our data were obtained for $\dot{\gamma}=10$ s$^{-1}$ whereas the data of Petrich et al. (2000) were measured at 0.5 s$^{-1}$. The fitting parameters are reported in Table II. $\sigma_{12}$ represents the shear stress component of the matrix.</figcaption></figure>

$$ \eta(\gamma_R) = \frac{4}{\gamma_R^4} \int_0^{\gamma_R} \eta(\gamma) \gamma^3 d\gamma, \qquad (32) $$

$$ N_1(\gamma_R) - N_2(\gamma_R) = \frac{2}{\gamma_R^2} \int_0^{\gamma_R} [N_1(\gamma) - N_2(\gamma)] \gamma d\gamma, \qquad (33) $$

where $\gamma_R$ is the strain at the disk rim and $\gamma$ is the deformation at the radius $r$. Furthermore, the IQUA [Eq. (28)] approximation is employed for the interaction tensor. The three quantities $X_A$, $k$, and $q$ are chosen as fitting parameters. A simple optimization technique was applied to determine the suitable coefficients, which were fixed for the four fiber concentrations.

Our reduced steady-state viscosity and normal stress difference data obtained at a shear rate of 10 s$^{-1}$ are compared in Fig. 8 to the data of Petrich et al. (2000) for similar,
<table><caption>**TABLE II.** Fitting parameters used to predict the shear viscosity and normal stress difference data of Fig. 8 and the transient data of Fig. 9.</caption><thead><tr><th></th><th>This work</th><th>Petrich <i>et al.</i><br>(2000)</th><th>Physical<br>meaning</th></tr></thead><tbody><tr><td><i>X</i><sub>A</sub></td><td>17.116</td><td>1.750</td><td>Parallel drag coefficient of the fiber</td></tr><tr><td><i>k</i></td><td>0.112</td><td>0.320</td><td>Intensity of the friction force</td></tr><tr><td><i>q</i></td><td>0.220</td><td>0.240</td><td>Interaction coefficient</td></tr></tbody></table>

but longer glass fiber suspensions, at a shear rate of 0.5 s⁻¹. In Petrich *et al.* (2000) the secondary normal stress differences were assumed to be negligible. It is interesting to note that both sets of data exhibit the same trend and, as expected, the rheological properties are larger for the larger aspect ratio fibers. The viscosity and normal stress data describe linear dependences at low fiber concentrations as observed by Sepehr *et al.* (2004c) for similar suspensions, but exhibit quadratic trends as fiber contents increase. The largest fiber volume fraction that was used by Sepehr *et al.* (2004c) was $\phi=7.06\%$, as compared to $\phi=10.5\%$ in our work. The quadratic dependences for both the reduced viscosity and normal stress differences are well predicted by the model. The parameters used to fit the data are presented in Table II. The same parameters have been used to fit the viscosity and normal stress data as well as the transient data presented in Fig. 9. We note that the parameters used to fit the data of Petrich *et al.* (2000) are largely different from the values to fit our data, although both suspensions are quite similar. This could be due to important differences in the fiber-fiber interactions for the two systems [the fibers used by Petrich *et al.* (2000) have been treated whereas in our case the treatment has been removed] or that our model overpredicts the effect of the aspect ratio [quadratic dependence on *r* predicted by Eq. (26)]. Kitano and Kataoka (1981) have also reported a quadratic behavior with volume fraction for vinyron and nylon fibers in silicone oil.

Figure 9 reports for our suspensions the transient shear viscosity and normal stress differences as functions of deformation for an applied shear rate of 10 s⁻¹. The stress growth data were obtained in reverse flow (CCW direction) following a first stress growth in the CW direction, since they have been observed to be more reproducible [Sepehr *et al.* (2004c)]. The data and model predictions for the four composites are presented in Fig. 9 using a semilog scale to highlight the behavior at low strain. In Fig. 9(a), we can see that the magnitude and the width of the stress growth overshoot, which have been related to fiber tumbling [Laun (1984); Sepehr *et al.* (2004b)], are fairly well described by the model using the same parameters (Table II). In Fig. 9(b), the model predictions for the normal stress differences are shown to be in good qualitative agreement with the data. The model predicts correctly that under reverse flow the normal stress differences take initially negative values (nevertheless underestimated) before rising and depicting a positive overshoot. However, the amplitude of the overshoots is slightly underestimated. Obviously, better fits would be obtained if different values for the parameters were used for each suspension.

Another interesting point is to compare the values obtained for the interaction coefficient using the new model with values obtained by Phan-Thien *et al.* (2002) by calculating the Folgar–Tucker constant $C_I$. The last term in Eq. (22) is a diffusive term for which we can write the interaction coefficient $C_I^*$ as

$$C_I^* = f \phi \tilde{M} q. \qquad (34)$$
<figure><img src="image_9.png" /><figcaption>FIG. 9. Stress growth data for our suspensions at |γ| = 10 s⁻¹ compared to model predictions for reverse flow (CCW) using the ORW3 and IQUA closure approximations. The fitting parameters are reported in Table II (a) Transient viscosity and (b) transient normal stress differences as functions of strain.</figcaption></figure>

As Eq. (22) contains another interaction term (third term on the right-hand side), $C_I^*$ is not really comparable to the Folgar–Tucker constant. The numerical values of $C_I^*$ for our data are about an order of magnitude smaller than values calculated by Phan-Thien *et al.* (2002) for $C_I$, but the trend shows a similar increase with fiber volume fraction. In contrast to Phan-Thien *et al.* expression [Eq. (3)], which is only dependent on the particle aspect ratio $r$ and the fiber volume fraction $\phi$, $C_I^*$ is proportional to $\phi$ and to $f$, which is a scalar function of the interfiber spacing, as suggested by Ranganathan and Advani (1991). It implies that the proposed diffusion coefficient evolves with time as $f$ depends on fiber orientation. This affects considerably the fiber orientation in transient flows and is of major importance for the modeling of industrial processes.

Finally, with the trace of the second-order interaction tensor, it is possible to predict the average contact number per fiber as a function of the material deformation using Eq. (21). This is given in Fig. 10 for reverse flow (CCW). As expected, a maximum of fiber interactions occurs at a given deformation and then the average fiber contacts decrease to
<figure><img src="image_10.png" /><figcaption>FIG. 10. Predictions of the average lubricated contact number per fiber as a function of deformation in reverse flow (CCW) at a shear rate of 10 s⁻¹ (ORW3 and IQUA closure approximations).</figcaption></figure>

reach a steady value. Moreover, the effective lubricated contact number per fiber is similar to the one obtained by Ausias *et al.* (2006) from direct simulations. For their largest concentrated suspensions ($\phi=11.5\%$), their computed average contact number was close to 2.0, similarly to what we find for $\phi=10.5\%$. To justify the trends in Fig. 10, we have computed the average spacing between fibers ($h$) according to Dinh and Armstrong (1984). For PB25, it is found that $h=(nL)^{-1/2} \approx 44$ µm for totally aligned fibers and $h=(nL^2)^{-1} \approx 6$ µm for isotropic orientation. Consequently, more contacts are expected when fibers disorient, as predicted by the model. The small oscillations, that appear at large deformation, as already observed in Figs. 2–5, are attributed to the used of the quadratic closure approximation for the fourth-order interaction tensor.

Our first observations show that the physics from which the model is derived is of a great relevance. The proposed approximations seem to be quite reasonable but further efforts are necessary, mainly concerning the closure approximations. The orientation probability distribution function should be numerically solved in order to suggest better closure approximations dealing with the new fourth-order interaction tensor. Improvement of the present work can also be accomplished by modeling the interacting force with a nonlinear lubrication expression such as a power law instead of a linear one. Finally, the obtained diffusivity is isotropic and improvement could be achieved by using an anisotropic diffusivity.

## VI. CONCLUSIONS

A new mathematical formulation has been proposed to describe the rheological behavior of fiber suspensions in the semiconcentrated regime and can be considered as an extension of the Dinh and Armstrong (1984) model. It takes into account fiber-fiber interactions with linear lubrication forces. In addition to the orientation tensors, this development allowed to introduce interaction tensors. A new equation of change for the second-order orientation tensor was derived, and an additional extra stress term was developed to take into account the contribution induced by the interacting forces. Closure approximations were proposed to estimate the fourth-order interaction tensors and the most accurate one is based on a quadratic derivation.
In order to validate the model, a rheological study of suspensions with fiber volume fractions close to those encountered in industrial applications have been performed using a Newtonian polybutene matrix and glass fibers. Despite the simple assumptions postulated at the fiber scale, start-up flow results showed a fairly good agreement in terms of the predicted viscosity and normal stress differences. The measured reduced viscosity and normal stress differences were found to be quadratic functions of the fiber volume fraction, as predicted by the model. Finally, the new model requires only a small set of constitutive parameters (3), which are all directly linked to the physics and geometry at the fiber scale. In the present study, the same three parameters have been used to predict all the investigated fiber suspensions at different volume fractions, ranging from the dilute to the semiconcentrated regimes.

## ACKNOWLEDGMENTS

This work was supported by the France-Québec collaboration program and funded by NSERC (Natural Sciences and Engineering Council of Canada)—CIAM (Inter-American Collaboration in Materials) program. The authors are also thankful to one reviewer who pointed out to them that the initially proposed model was not objective.

## APPENDIX A: REDUCED-STRAIN CLOSURE MODEL

Recently, Tucker et al. (2007) and Phelps et al. (2008) proposed a new evolution equation for the second orientation tensor. The so-called reduced-strain closure (RSC) model slows down the orientation kinetics and is given for an infinite fiber ratio by

$$ \frac{Da_{ij}^{\text{RSC}}}{Dt} = -\frac{1}{2}(\omega_{ik}a_{kj} - a_{ik}\omega_{kj}) + \frac{1}{2}\{\dot{\gamma}_{ik}a_{kj} + a_{ik}\dot{\gamma}_{kj} - 2[a_{ijkl} + (1-\kappa)(L_{ijkl} - M_{ijmn}a_{nmkl})]\dot{\gamma}_{kl}\} + 2\kappa C_i |\dot{\gamma}|(\delta_{ij} - 3a_{ij}), \quad (A1) $$

where **L** and **M** are fourth-order tensors, which are calculated from eigenvalues and eigenvectors of the second-order orientation tensor.

By factorizing the deformation rate tensor $\dot{\gamma}$ in the third member on the right-hand side of Eq. (22), the evolution equation for the second order orientation tensor can be rewritten as

$$ \frac{Da_{ij}}{Dt} = -\frac{1}{2}(\omega_{ik}a_{kj} - a_{ik}\omega_{kj}) + \frac{1}{2}\left(\dot{\gamma}_{ik}a_{kj} + a_{ik}\dot{\gamma}_{kj} - 2\left\{a_{ijkl} + \phi \tilde{M} \left[\frac{1}{4}(b_{ik}\delta_{jl} + b_{jl}\delta_{ik} + b_{il}\delta_{jk} + b_{jk}\delta_{il}) - b_{ijkl}\right]\right\}\dot{\gamma}_{kl}\right) + 2f\phi \tilde{M} q |\dot{\gamma}|(\delta_{ij} - 3a_{ij}). \quad (A2) $$

This shows that our equation for the time derivative of $a_2$ is structured as a reduced-strain closure model.

## APPENDIX B: FIBER MASS CENTER VELOCITY

The global hydrodynamic force $F_H$ applied along the fiber length in the dilute regime is expressed as
$$F_H = A \cdot (v^\infty - \dot{r}_C^\alpha), \quad (B1)$$

where $\dot{r}_C^\alpha$ is the mass center velocity of the test particle. A is a second order hydrodynamic resistance tensor, which contains fluid drag coefficients and is expressed by $A = \eta_0 L[Y_A \delta+(X_A-Y_A)\mathbf{p}^\alpha\mathbf{p}^\alpha]$ in the R coordinate systems (Fig. 1). Coefficients $X_A$ and $Y_A$ describe the drag exerted on the fiber from the axial and transverse fluid motion, respectively. For dilute suspensions of fibers in a Newtonian fluid, the slender body theory results in $Y_A=2X_A=2$ [Dinh and Armstrong (1984)]. A force balance for negligible fiber inertia gives $\dot{r}_C^\alpha = \kappa \cdot \mathbf{r}^\alpha$ as derived by Dinh and Armstrong (1984), and it shows that the fiber centroid moves affinely with the effective medium in the absence of any other forces besides those imparted by the suspending fluid.

Furthermore, in the case where contacts between fibers are taken into account, the total force on the test fiber is now the sum of the global hydrodynamic and interaction force. Again, a force balance for negligible fiber inertia gives the motion of the center mass

$$\dot{\mathbf{r}}_C^\alpha = \kappa \cdot \mathbf{r}_C^\alpha + 2k \eta_0 n D^2 L^2 \mathbf{A}^{-1} \cdot (\dot{\mathbf{r}}_C^\beta - \dot{\mathbf{r}}_C) \int_{\mathbf{p}^\beta} |\mathbf{p}^\alpha \times \mathbf{p}^\beta| \psi_{\mathbf{p}^\beta} d\mathbf{p}^\beta, \quad (B2)$$

where $\dot{\mathbf{r}}_C^\beta$ is the mass center velocity of the neighboring fiber. Consequently, the center mass deforms non-affinely if the interactions are present. The Dinh and Armstrong (1984) expression is recovered when the number of fibers per unit volume (n) is small.

However, if we suppose that $\dot{\mathbf{r}}_C^\beta$ is close to $\dot{\mathbf{r}}_C^\alpha$ plus a perturbation velocity, as done for $\dot{\mathbf{p}}^\alpha$ in Appendix C, the next result is obtained after simplification

$$\dot{\mathbf{r}}_C^\alpha = \kappa \cdot \mathbf{r}_C^\alpha - 2knD^2L \frac{D_t}{Y_A \psi_{\mathbf{p}^\alpha}} \mathbf{A}^{-1} \cdot \frac{\partial \psi_{\mathbf{p}^\alpha}}{\partial \mathbf{r}_C^\alpha} \int_{\mathbf{p}^\beta} |\mathbf{p}^\alpha \times \mathbf{p}^\beta| \psi_{\mathbf{p}^\beta} d\mathbf{p}^\beta, \quad (B3)$$

where the parameter $D_t$ is expressed in m/s. According to hypothesis (3), it implies that

$$\dot{\mathbf{r}}_C^\alpha = \kappa \cdot \mathbf{r}_C^\alpha. \quad (B4)$$

This result clearly shows that the mass center of the fibers moves affinely with the bulk fluid.

## APPENDIX C: EVALUATION OF INTERACTION TENSORS WITH MAIER-SAUPE APPROXIMATION

The scalar quantity can be extended as follows:

$$\begin{align}
|\mathbf{p}^\alpha \times \mathbf{p}^\beta| &\Rightarrow |\mathbf{p}^\alpha \times \mathbf{p}^\beta|^2 = (\mathbf{p}^\alpha \times \mathbf{p}^\beta)^2 = (\mathbf{p}^\alpha \times \mathbf{p}^\beta)_i (\mathbf{p}^\alpha \times \mathbf{p}^\beta)_i = (\varepsilon_{ijk} p_j^\alpha p_k^\beta)(\varepsilon_{imn} p_m^\alpha p_n^\beta) \nonumber \\
&= (\varepsilon_{ijk} \varepsilon_{imn})(p_j^\alpha p_k^\beta p_m^\alpha p_n^\beta) = (\delta_{jm} \delta_{kn} - \delta_{jn} \delta_{km})(p_j^\alpha p_k^\beta p_m^\alpha p_n^\beta) \nonumber \\
&= \delta_{jm} \delta_{kn} p_j^\alpha p_k^\beta p_m^\alpha p_n^\beta - \delta_{jn} \delta_{km} p_j^\alpha p_k^\beta p_m^\alpha p_n^\beta = p_j^\alpha p_k^\beta - p_j^\alpha p_k^\beta - p_j^\alpha p_k^\beta = 1 - \mathbf{p}^\alpha \mathbf{p}^\alpha : \mathbf{p}^\beta \mathbf{p}^\beta. \tag{C1}
\end{align}$$

Moreover, we introduce a numerical prefactor of $3\pi/8$ to reach the exact value of $\pi/4$ for the isotropic average value. After replacing expression (C1) into Eq. (16), we get
$$
\begin{align*}
\mathbf{b}_2 &= \int_{\mathbf{p}^\alpha} \int_{\mathbf{p}^\beta} \mathbf{p}^\alpha \mathbf{p}^\alpha |\mathbf{p}^\alpha \times \mathbf{p}^\beta| \psi_{\mathbf{p}^\alpha} \psi_{\mathbf{p}^\beta} d\mathbf{p}^\alpha d\mathbf{p}^\beta \\
&= \frac{3\pi}{8} \int_{\mathbf{p}^\alpha} \int_{\mathbf{p}^\beta} \mathbf{p}^\alpha \mathbf{p}^\alpha (1 - \mathbf{p}^\alpha \mathbf{p}^\alpha : \mathbf{p}^\beta \mathbf{p}^\beta) \psi_{\mathbf{p}^\alpha} \psi_{\mathbf{p}^\beta} d\mathbf{p}^\alpha d\mathbf{p}^\beta \\
&= \frac{3\pi}{8} \int_{\mathbf{p}^\alpha} \int_{\mathbf{p}^\beta} (\mathbf{p}^\alpha \mathbf{p}^\alpha - \mathbf{p}^\alpha \mathbf{p}^\alpha \mathbf{p}^\alpha : \mathbf{p}^\beta \mathbf{p}^\beta) \psi_{\mathbf{p}^\alpha} \psi_{\mathbf{p}^\beta} d\mathbf{p}^\alpha d\mathbf{p}^\beta \\
&= \frac{3\pi}{8} \int_{\mathbf{p}^\alpha} (\mathbf{p}^\alpha \mathbf{p}^\alpha - \mathbf{p}^\alpha \mathbf{p}^\alpha \mathbf{p}^\alpha : \mathbf{a}_2) \psi_{\mathbf{p}^\alpha} d\mathbf{p}^\alpha = \frac{3\pi}{8} (\mathbf{a}_2 - \mathbf{a}_4 : \mathbf{a}_2). \tag{C2}
\end{align*}
$$

References

Advani, S. G., and C. L. Tucker, "The use of tensors to describe and predict fiber orientation in short fiber composites," J. Rheol. 31, 751-784 (1987).

Ausias, G., J. F. Agassant, and M. Vincent, "Flow and fiber orientation calculations in reinforced thermoplastic extruded tubes," Int. Polym. Process. 9, 51-59 (1994).

Ausias, G., X. J. Fan, and R. I. Tanner, "Direct simulation for concentrated fibre suspensions in transient and steady state shear flows," J. Non-Newtonian Fluid Mech. 135, 46-57 (2006).

Batchelor, G. K., "The stress generated in a non-dilute suspension of elongated particles by pure straining motion," J. Fluid Mech. 46, 813-829 (1971).

Batchelor, G. K., "The stress system in a suspension of force-free particles," J. Fluid Mech. 41, 545-570 (1970).

Bay, R. S., "Fiber orientation in injection-molded composites: A comparison of theory and experiment," Ph.D. dissertation, University of Illinois at Urbana-Champaign, IL, 1991.

Bird, R. B., C. F. Curtiss, R. C. Armstrong, and O. Hassager, *Dynamics of polymeric liquids. Volume 2, Kinetic theory* (Wiley, New York, 1987).

Carreau, P. J., D. De Kee, and R. P. Chhabra, *Rheology of Polymeric Systems: Principles and Applications* (Hanser, Munich, 1997).

Chaouche, M., and D. L. Koch, "Rheology of non-Brownian rigid fiber suspensions with adhesive contacts," J. Rheol. 45, 369-382 (2001).

Chung, D. H., and T. H. Kwon, "Fiber orientation in the processing of polymer composites," Korea-Aust. Rheol. J. 14, 175-188 (2002).

Chung, D. H., and T. H. Kwon, "Improved model of orthotropic closure approximation for flow induced fiber orientation," Polym. Compos. 22, 636-649 (2001).

Claeys, I. L., and J. F. Brady, "Suspensions of prolate spheroids in Stokes flow. Part 1. Dynamics of finite number of particles in an unbounded fluid," J. Fluid Mech. 251, 411-442 (1993).

Dhont, J. K. G., and W. J. Briels, "Viscoelasticity of suspensions of long, rigid rods," Colloids Surf., A 213, 131-156 (2003).

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

Djalili-Moghaddam, M., and S. Toll, "A model for short-range interactions in fibre suspensions," J. Non-Newtonian Fluid Mech. 132, 73-83 (2005).

Doi, M., and S. F. Edwards, "Dynamics of rod-like macromolecules in concentrated solution. Part 1," J. Chem. Soc., Faraday Trans. 2 74, 560-570 (1978).

Folgar, F. P., and C. L. Tucker, "Orientation behavior of fibers in concentrated suspensions," J. Reinf. Plast. Compos. 3, 99-119 (1984).

Gibson, A. G., and S. Toll, "Mechanics of the squeeze flow of planar fibre suspensions," J. Non-Newtonian
Fluid Mech. **82**, 1–24 (1999).

Giesekus, H., "Elasto-viskose flüssigkeiten, für die in stationären schichtströmungen sämtliche normalspannungskomponenten verschieden gross sind," Rheol. Acta **2**, 50–62 (1962).

Grmela, M., and B. Z. Dlugogorski, "Hamiltonian modelling of liquid crystal polymers and blends," in *Rheology and Processing of Liquid Crystal Polymers* (Chapman & Hall, London, 1996).

Hand, G. L., "A theory of dilute suspensions," Arch. Ration. Mech. Anal. **7**, 81–86 (1961).

Jeffery, G. B., "The motion of ellipsoidal particles immersed in a viscous fluid," Proc. R. Soc. London, Ser. A **102**, 161–179 (1922).

Joung, C. G., N. Phan-Thien, and X. J. Fan, "Direct simulation of flexible fibers," J. Non-Newtonian Fluid Mech. **99**, 1–36 (2001).

Kim, S., and S. J. Karrila, *Microdynamics: Principles and Selected Applications* (Butterworth-Heinemann, Boston, 1991).

Kitano, T., and T. Kataoka, "The rheology of suspensions of vinyloxycarbonyl fibers in polymer liquids. I. Suspensions in silicone oil," Rheol. Acta **20**, 390–402 (1981).

Laun, H. M., "Orientation effects and rheology of short glass fiber-reinforced thermoplastics," Colloid Polym. Sci. **262**, 257–269 (1984).

Le Corre S., P. Dumont, L. Orgeas, and D. Favier, "Rheology of highly concentrated planar fiber suspensions," J. Rheol. **49**, 1029–1058 (2005).

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

Maier, W., and A. Saupe, "Eine einfache molekular-statistische Theorie des nematischen kristallinflüssigen Zustandes," Z. Naturforsch. A **14A**, 882–889 (1959).

Onsager, L., "The effects of shape on the interaction of colloidal particles," Ann. N.Y. Acad. Sci. **51**, 627–659 (1949).

Peterlin, A., "Excluded volume effect on light scattering of the coiled linear macromolecule," J. Chem. Phys. **23**, 2464–2465 (1955).

Petrich, M. P., D. L. Koch, and C. Cohen, "An experimental determination of the stress-microstructure relationship in semi-concentrated fiber suspensions," J. Non-Newtonian Fluid Mech. **95**, 101–133 (2000).

Phan-Thien, N., X. J. Fan, R. I. Tanner, and R. Zheng, "Folgar-Tucker constant for a fibre suspension in a Newtonian fluid," J. Non-Newtonian Fluid Mech. **103**, 251–260 (2002).

Phelps, J. H., J. Wang, and C. L. Tucker, "New model for fiber orientation in injection-molded composites," in 24th Annual Meeting of the Polymer Processing Society, Salerno, Italy, 2008.

Ranganathan, S., and S. G. Advani, "Fiber-fiber interactions in homogeneous flows of nondilute suspensions," J. Rheol. **35**, 1499–1522 (1991).

Sandstrom, C. R., "Interactions and orientation in concentrated suspensions of rigid rods: Theory and experiment," Ph.D. thesis, University of Illinois at Urbana-Champaign, IL, 1993.

Sandstrom, C. R., and C. L. Tucker, "A theory for concentrated fiber suspensions with strong fiber-fiber interactions," Makromol. Chem., Macromol. Symp. **68**, 291–300 (1993).

Schmid, C. F., L. H. Switzer, and D. J. Klingenberg, "Simulations of fiber flocculation: Effects of fiber properties and interfiber friction," J. Rheol. **44**, 781–809 (2000).

Sepehr, M., G. Ausias, and P. J. Carreau, "Rheological behavior in transient shear flows for concentrated short fiber filled polymer," Mec. Ind. **5**, 419–428 (2004a).

Sepehr, M., G. Ausias, and P. J. Carreau, "Rheological properties of short fiber filled polypropylene in transient shear flow," J. Non-Newtonian Fluid Mech. **123**, 19–32 (2004b).

Sepehr, M., P. J. Carreau, M. Moan, and G. Ausias, "Rheological properties of short fiber model suspensions," J. Rheol. **48**, 1023–1048 (2004c).

Servais, C., A. Luciani, and J. A. E. Manson, "Fiber-fiber interaction in concentrated suspensions: Dispersed fiber bundles," J. Rheol. **43**, 1005–1018 (1999a).

Servais, C., J. A. E. Manson, and S. Toll, "Fiber-fiber interaction in concentrated suspensions: Disperse fibers," J. Rheol. **43**, 991–1004 (1999b).

Shaqfeh, E. S. G., and G. H. Fredrickson, "The hydrodynamic stress in a suspension of rods," Phys. Fluids A **2**, 7–24 (1990).
Sundararajakumar, R. R., and D. L. Koch, "Structure and properties of sheared fiber suspensions with mechanical contacts," J. Non-Newtonian Fluid Mech. **73**, 205-239 (1997).

Toll, S., "Packing mechanics of fiber reinforcements," Polym. Eng. Sci. **38**, 1337-1350 (1998).

Toll, S., and J. A. E. Manson, "Dynamics of a planar concentrated fiber suspension with nonhydrodynamic interaction," J. Rheol. **38**, 985-997 (1994).

Tucker, C. L., J. Wang, and J. F. O'Gara, "Method and article of manufacture for determining a rate of change of orientation of a plurality of fibers disposed in a fluid," U.S. Patent No. 7,266,469 B1 (2007).

Yamamoto, S., and T. Matsuoka, "Dynamic simulation of microstructure and rheology of fiber suspensions," Polym. Eng. Sci. **36**, 2396-2403 (1996).