2.1 Cell model

The test rod is coaxially embedded into a cylindrical fluid cell of radius $h$ and of the same length as the particle (figure 1). This cell model assumes that the effect of one fibre on its neighbours is approximated by an equivalent cylindrical boundary around the test rod, and therefore simplifies the problem to a single-particle theory. In order to determine the velocity field near a slender particle, a dimensional analysis is carried out using a cylindrical coordinate system for convenience. First, the fluid velocity field is considered in the vicinity ( $|z|<L/2$ and $r\ll L/2$ ) of a particle in the assumed equilibrium orientation, as shown schematically in figure 1. Provided that the cross-sectional form of the particle is circular and does not vary with the axial $z$ -position, it is useful to adopt a dimensional scaling for the velocities and the velocity gradients. Thus, the dimensionless quantities are distinguished by asterisks such as

(2.1a,b ) $$\begin{eqnarray}r=Rr^{\ast },\quad {\it\theta}={\it\theta}^{\ast },\quad z=\frac{L}{2}z^{\ast }\quad \text{and}\quad u_{r}=R|\dot{{\it\gamma}}|u_{r}^{\ast },\quad u_{{\it\theta}}=R|\dot{{\it\gamma}}|u_{{\it\theta}}^{\ast },\quad u_{z}=\frac{L}{2}|\dot{{\it\gamma}}|u_{z}^{\ast }.\end{eqnarray}$$

We are concerned with the limiting forms of the dimensionless velocities $u_{r}^{\ast }$ , $u_{{\it\theta}}^{\ast }$ and $u_{z}^{\ast }$ and their $r^{\ast }$ , ${\it\theta}^{\ast }$ , $z^{\ast }$ -derivatives for $a_{r}\rightarrow \infty$ in an inner region near the body, where $r^{\ast }$ , ${\it\theta}^{\ast }$ and $z^{\ast }$ are all $O(1)$ . Since the gradient operator is given by

(2.2) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}={\bf\delta}_{r}\frac{\partial }{\partial r}+{\bf\delta}_{{\it\theta}}\frac{1}{r}\frac{\partial }{\partial {\it\theta}}+{\bf\delta}_{z}\frac{\partial }{\partial z}=\frac{1}{R}\left(\boldsymbol{{\rm\nabla}}^{(0)}+\frac{1}{a_{r}}\boldsymbol{{\rm\nabla}}^{(1)}\right),\end{eqnarray}$$

where

(2.3a,b ) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}^{(0)}={\bf\delta}_{r}\frac{\partial }{\partial r^{\ast }}+{\bf\delta}_{{\it\theta}}\frac{1}{r^{\ast }}\frac{\partial }{\partial {\it\theta}^{\ast }}\quad \text{and}\quad \boldsymbol{{\rm\nabla}}^{(1)}={\bf\delta}_{z}\frac{\partial }{\partial z^{\ast }},\end{eqnarray}$$

one has formally for the velocity gradient

(2.4) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{u}=a_{r}|\dot{{\it\gamma}}|[\boldsymbol{{\rm\nabla}}^{(0)}({\bf\delta}_{z}u_{z}^{\ast })]=a_{r}|\dot{{\it\gamma}}|\left[{\bf\delta}_{r}{\bf\delta}_{z}\frac{\partial u_{z}^{\ast }}{\partial r^{\ast }}+O(u/a_{r})\right],\end{eqnarray}$$

where the ${\it\theta}^{\ast }$ -derivatives have been omitted due to the problem symmetry. Here, we denote by $O(u/a_{r})$ the terms which are $O(1/a_{r})$ relative to the leading terms in (2.4). The leading terms themselves suggest that a quasi-steady and shear-dominated flow occurs in the near field, which implies the validity of a viscometric flow representation for the fluid rheology (Batchelor Reference Batchelor1971; Goddard Reference Goddard1976a ,Reference Goddard b ). Therefore, this cell model provides a useful approximation in the particle near field, with constant velocity gradients along the rod.

Figure 1. Geometry of the cell surrounding the test rod.

For this axisymmetric problem, the axial velocity $u_{z}$ satisfies the following equation of motion in terms of stress in cylindrical coordinates

(2.5) $$\begin{eqnarray}\frac{1}{r}\frac{\partial }{\partial r}(r{\it\tau}_{rz})=0,\end{eqnarray}$$

where ${\it\tau}_{rz}$ represents the shear stress. First, equation (2.5) can be integrated directly to give the familiar elementary form for the shear stress ${\it\tau}_{rz}={\it\tau}_{R}R/r$ , where ${\it\tau}_{R}$ is the unknown stress at the rod surface (a no-slip condition is considered at the particle surface). Secondly, it is assumed that the rheological relationship can be inverted to give the shear rate, $\dot{{\it\gamma}}$ , as a function of the shear stress. Hence, following a change of variable (i.e. $r=R{\it\tau}_{R}/{\it\tau}_{rz}$ ), the axial velocity $u_{z}$ can be expressed as

(2.6) $$\begin{eqnarray}u_{z}=\int _{R}^{r}\dot{{\it\gamma}}\,\text{d}r={\it\tau}_{R}R\int _{{\it\tau}_{R}R/r}^{{\it\tau}_{R}}\frac{\dot{{\it\gamma}}}{{\it\tau}_{rz}^{2}}\,\text{d}{\it\tau}_{rz}.\end{eqnarray}$$

Once the shear rheology of the fluid is known, this result provides a relationship between the axial velocity distribution and the rod surface stress ${\it\tau}_{R}$ . At the same time, it becomes evident that an outer boundary condition is required to determine ${\it\tau}_{R}$ . At $r=h$ , the relative fluid velocity field is assumed to be undisturbed by the particle. For an inertialess particle, a force balance shows that its centre of mass translates affinely with the bulk flow. Hence, the fluid velocity at distance $z\boldsymbol{p}$ from the centre of mass of the rod is $\boldsymbol{u}^{fluid}=z{\bf\kappa}\boldsymbol{\cdot }\boldsymbol{p}$ , where ${\bf\kappa}$ is the transpose of the velocity gradient tensor (Bird et al. Reference Bird, Armstrong and Hassager1987a ), whereas the rod velocity at the same location is $\boldsymbol{u}^{rod}=z\dot{\boldsymbol{p}}$ . The relative velocity is $\boldsymbol{u}=\boldsymbol{u}^{fluid}-\boldsymbol{u}^{rod}$ , and hence its component along the $z$ -direction, is found to be $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{p}=z{\bf\kappa}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}$ , which gives the condition $u_{z}(h)=z\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}/2$ (when expressed in the local coordinate system). $z$ is the arc length measured along the rod axis of direction $\boldsymbol{p}$ with $z=0$ at the centre of the rod (figure 1). With this result in mind, equation (2.6) yields

(2.7) $$\begin{eqnarray}z\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}/2={\it\tau}_{R}R\int _{{\it\tau}_{R}R/h}^{{\it\tau}_{R}}\frac{\dot{{\it\gamma}}}{{\it\tau}_{rz}^{2}}\,\text{d}{\it\tau}_{rz}.\end{eqnarray}$$

Note that $\dot{{\bf\gamma}}$ is the unperturbed deformation rate tensor taken at the outer boundary of the cell and $\dot{{\it\gamma}}=\partial u_{z}/\partial r$ is the disturbed shear rate close to the rod. Equation (2.7) is a key formula used throughout this article in order to take into account the contribution of particles suspended in nonlinear matrices. An elementary force balance on the particle leads to the force per unit length $f=2{\rm\pi}R{\it\tau}_{R}$ . This axial tensile force is parallel to the unit vector $\boldsymbol{p}$ , as Batchelor (Reference Batchelor1971) argued that the components normal to the particle make no contribution in the case of a rod on which no external force or couple acts.

The average force carried by the rod is then substituted into the Kramer expression (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987b ) to derive the rod stress tensor, in which the volume fraction of particles, ${\it\phi}$ , takes into account the contribution of each particle contained in a given volume. Hence, the extra stress due to the particle contribution is

(2.8) $$\begin{eqnarray}{\bf\tau}^{p}=\frac{2{\it\phi}}{{\rm\pi}R^{2}L}\left\langle \boldsymbol{p}\boldsymbol{p}\int _{z=0}^{z=L/2}zf\,\text{d}z\right\rangle ,\end{eqnarray}$$

where the angular brackets denote the ensemble average with respect to the distribution function of $\boldsymbol{p}$ . Equation (2.8) shows that the stress can be directly computed once the drag force is known.

2.2 Results for Newtonian fluids

If a Newtonian medium of viscosity ${\it\eta}_{0}$ is considered, the usual rheological constitutive equation expressed in simple shear flow is $\dot{{\it\gamma}}={\it\tau}_{rz}/{\it\eta}_{0}$ and (2.7) becomes

(2.9) $$\begin{eqnarray}f=\frac{{\rm\pi}{\it\eta}_{0}}{\log (h/R)}z\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p},\end{eqnarray}$$

which is the well-known parallel drag force for a rod suspended into a Newtonian fluid. Combining (2.8) and (2.9) leads to the following particle contribution to the stress tensor

(2.10) $$\begin{eqnarray}{\bf\tau}^{p}={\it\eta}_{0}{\it\phi}\frac{a_{r}^{2}}{3\log (h/R)}\dot{{\bf\gamma}}\boldsymbol{ : }\langle \boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\rangle .\end{eqnarray}$$

In order to derive the evolution equation for the second-order tensor $\langle \boldsymbol{p}\boldsymbol{p}\rangle$ , Giesekus (Reference Giesekus1962) proposed an ingenious method to determine the stress tensor (which is frequently used in the kinetic theory of polymer solutions and melts (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987b )) and is referred to as the Giesekus stress tensor. This result can be used to derive the extra stress due to the particle contribution and is given by

(2.11) $$\begin{eqnarray}{\bf\tau}^{p}=-\frac{n{\it\zeta}}{4}\frac{\mathscr{D}\langle \boldsymbol{p}\boldsymbol{p}\rangle }{\mathscr{D}t},\end{eqnarray}$$

where $n$ represents the number of particles per unit volume. In order to define this derivative, we remark that the Giesekus formulation involves a convected derivative of the contravariant components of a second-order tensor, ${\it\bf\Delta}$ , defined by

(2.12) $$\begin{eqnarray}\frac{\mathscr{D}{\it\bf\Delta}}{\mathscr{D}t}=\frac{\text{D}{\it\bf\Delta}}{\text{D}t}+\frac{1}{2}({\bf\omega}\boldsymbol{\cdot }{\it\bf\Delta}-{\it\bf\Delta}\boldsymbol{\cdot }{\bf\omega})-\frac{1}{2}(\dot{{\bf\gamma}}\boldsymbol{\cdot }{\it\bf\Delta}+{\it\bf\Delta}\boldsymbol{\cdot }\dot{{\bf\gamma}}).\end{eqnarray}$$

A comparison of (2.10) and (2.11) with the help of (2.12) leads to

(2.13) $$\begin{eqnarray}\frac{\text{D}\langle \boldsymbol{p}\boldsymbol{p}\rangle }{\text{D}t}=-\frac{1}{2}({\bf\omega}\boldsymbol{\cdot }\langle \boldsymbol{p}\boldsymbol{p}\rangle -\langle \boldsymbol{p}\boldsymbol{p}\rangle \boldsymbol{\cdot }{\bf\omega})+\frac{1}{2}(\dot{{\bf\gamma}}\boldsymbol{\cdot }\langle \boldsymbol{p}\boldsymbol{p}\rangle +\langle \boldsymbol{p}\boldsymbol{p}\rangle \boldsymbol{\cdot }\dot{{\bf\gamma}}-2\dot{{\bf\gamma}}\boldsymbol{ : }\langle \boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\rangle ),\end{eqnarray}$$

and the drag coefficient is found to be

(2.14) $$\begin{eqnarray}{\it\zeta}=\frac{{\it\eta}_{0}{\rm\pi}L^{3}}{3\log (h/R)}.\end{eqnarray}$$

These last results are in agreement with those obtained by Dinh & Armstrong (Reference Dinh and Armstrong1984) for semi-concentrated fibre suspensions in Newtonian fluids. Hence, a rheological constitutive equation for slender rods in Newtonian fluids can be expressed as follows

(2.15) $$\begin{eqnarray}{\bf\sigma}=-P{\bf\delta}+{\it\eta}_{0}\dot{{\bf\gamma}}+{\it\eta}_{0}{\it\phi}\frac{a_{r}^{2}}{3\log (h/R)}\dot{{\bf\gamma}}\boldsymbol{ : }\langle \boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\rangle .\end{eqnarray}$$

3.1 Constitutive equation for suspensions in power-law fluids

For power-law fluids, the shear-rate versus shear-stress rheological relation can be written as $\dot{{\it\gamma}}={\it\tau}_{rz}/K|{\it\tau}_{rz}/K|^{(1-m)/m}$ . The latter expression is substituting in (2.7) to give the following drag force

(3.1) $$\begin{eqnarray}f=\frac{{\rm\pi}K}{(2R)^{m-1}}\left[\frac{1-m}{m[1-(R/h)^{(1-m)/m}]}\right]^{m}z^{m}|\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}.\end{eqnarray}$$

It is interesting to note the nonlinear dependency of the drag force on the magnitude of the strain-rate tensor and the particle orientation through the dyadic product of $\boldsymbol{p}$ . This expression is then inserted in (2.8) and an analytical integration along the arc length can be performed to result in

(3.2) $$\begin{eqnarray}{\bf\tau}^{p}=K{\it\phi}\frac{a_{r}^{m+1}}{2^{m-1}(m+2)}\left[\frac{1-m}{m[1-(R/h)^{(1-m)/m}]}\right]^{m}\dot{{\bf\gamma}}\boldsymbol{ : }\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\rangle ,\end{eqnarray}$$

which is exactly the stress contribution for rod suspensions in power-law fluids obtained by Souloumiac & Vincent (Reference Souloumiac and Vincent1998) and Gibson & Toll (Reference Gibson and Toll1999). In the special case of $m=1$ , the solution reduces to (2.10) with $K={\it\eta}_{0}$ . Examination of (3.2) reveals that the slope of log ${\it\tau}_{12}^{p}$ versus $\log \dot{{\it\gamma}}$ , where ${\it\tau}_{12}^{p}$ is the particle shear stress, is given by the exponent of $m$ , which is the same slope as the suspending fluid: the model predicts that the presence of rods does not alter the slope with which the suspension shear thins. A rheological constitutive equation for axisymmetric particles in power-law fluids can then be represented as follows

(3.3) $$\begin{eqnarray}{\bf\sigma}=-P{\bf\delta}+K|\dot{{\it\gamma}}|^{m-1}\dot{{\bf\gamma}}+K{\it\phi}\frac{a_{r}^{m+1}}{2^{m-1}(m+2)}\left[\frac{1-m}{m[1-(R/h)^{(1-m)/m}]}\right]^{m}\dot{{\bf\gamma}}\boldsymbol{ : }\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\rangle .\end{eqnarray}$$

To derive the time evolution equation for the rod microstructure, the Giesekus expression is used again and yields to the following expression

(3.4) $$\begin{eqnarray}\displaystyle \frac{D\langle |\dot{{\bf\gamma}}\mathbf{:}\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle }{\text{D}t} & = & \displaystyle -\frac{1}{2}({\bf\omega}\boldsymbol{\cdot }\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle -\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle \boldsymbol{\cdot }{\bf\omega})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{2} (\dot{{\bf\gamma}}\boldsymbol{\cdot }\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle +\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle \boldsymbol{\cdot }\dot{{\bf\gamma}}\nonumber\\ \displaystyle & & \displaystyle -\,2\dot{{\bf\gamma}}\boldsymbol{ : }\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\rangle ),\end{eqnarray}$$

where $\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle$ is a symmetric second-order tensor. In order to gain insight into the orientation evolution of a rod, the time evolution equation for the probability distribution function, assuming the existence of a diffusive contribution as suggested by Folgar & Tucker (Reference Folgar and Tucker1984), is given by

(3.5) $$\begin{eqnarray}\frac{\text{D}{\it\psi}}{\text{D}t}=-\frac{\partial }{\partial \boldsymbol{p}}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}{\it\psi})+C_{I}|\dot{{\it\gamma}}|\frac{\partial ^{2}{\it\psi}}{\partial \boldsymbol{p}^{2}},\end{eqnarray}$$

where $\partial /\partial \boldsymbol{p}$ is the differential operator on the surface of a unit sphere. $\dot{\boldsymbol{p}}$ remains to be expressed considering that the suspending fluid follows a power-law behaviour. Hence, it is assumed at this point that $\dot{\boldsymbol{p}}$ is given by the Jeffery equation (1.3) for a slender body (i.e. ${\it\lambda}=1$ ) with no diffusion (i.e. $C_{I}=0$ ). Although Jeffery (Reference Jeffery1922) solved the velocity field around a single rigid ellipsoidal particle immersed in a Newtonian fluid with negligible inertia, this assumption is justified as for a single material point in simple shear flow, the Newtonian viscosity of the surrounding fluid can be replaced by a local viscosity predicted by the power-law model (still a viscous fluid with no time dependence and no normal stress differences). This assumption appears to be all the more valid since the viscosity term does not appear in the Jeffery solution. In our cell model, the velocity gradient is also assumed constant along the length of the particle. Moreover, a force balance using (3.1) leads to the rod centroid moving affinely with regards to the effective medium in the absence of any other forces besides those imparted by the suspending fluid, as also suggested by Jeffery (Reference Jeffery1922). Thus the Jeffery equation is substituted into (3.5) and the convective part for the time evolution of $\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle$ is found to be (3.4) (details on the derivation are given in appendix A). We were unfortunately unable to straightforwardly express the time evolution for the diffusive term.

Dividing (3.4) by $\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\rangle$ , which corresponds to the trace of $\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle$ , equation (3.4) results in

(3.6) $$\begin{eqnarray}\frac{\text{D}\unicode[STIX]{x1D656}_{2}^{(m)}}{\text{D}t}=-\frac{1}{2}({\bf\omega}\boldsymbol{\cdot }\unicode[STIX]{x1D656}_{2}^{(m)}-\unicode[STIX]{x1D656}_{2}^{(m)}\boldsymbol{\cdot }{\bf\omega})+\frac{1}{2}(\dot{{\bf\gamma}}\boldsymbol{\cdot }\unicode[STIX]{x1D656}_{2}^{(m)}+\unicode[STIX]{x1D656}_{2}^{(m)}\boldsymbol{\cdot }\dot{{\bf\gamma}}-2\dot{{\bf\gamma}}\boldsymbol{ : }\unicode[STIX]{x1D656}_{4}^{(m)}),\end{eqnarray}$$

where $\unicode[STIX]{x1D656}_{2}^{(m)}=\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\rangle /\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\rangle$ and $\unicode[STIX]{x1D656}_{4}^{(m)}=\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\boldsymbol{p}\rangle /\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}\rangle$ have been used. Similarly to the orientation tensors (Advani & Tucker Reference Advani and Tucker1987), $\unicode[STIX]{x1D656}_{2}^{(m)}$ and $\unicode[STIX]{x1D656}_{4}^{(m)}$ describe an equivalent orientation state for the rods, for which the drift due to the scalar potential $|\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}$ is considered (note that the trace of $\unicode[STIX]{x1D656}_{2}^{(m)}$ always equals 1).

Equation (3.6) is structured as (1.4) which is expressed for a Newtonian matrix with infinite particle aspect ratio ( ${\it\lambda}=1$ ) and no diffusion ( $C_{I}=0$ ). Therefore, these compact and efficient tensors are also called orientation tensors. The components of the second-order orientation tensor have some simple physical interpretations. If all its diagonal components are $1/3$ , the orientation tensor indicates a three-dimensional isotropic orientation state for the rods. When two diagonal components are equal to $1/2$ , the tensor implies a two-dimensional random, planar bi-axial or planar orientation. Finally, a perfect uniaxial alignment results in one diagonal component value of 1.

The previous governing equations (3.3) and (3.4) are not satisfactory to completely close the problem: one needs some approximations such as closure relations. The objective of this work is not to focus on the development of such relations and the accompanying issue of questionable accuracy. Therefore, equation (3.6) will no longer be used in the rest of the paper and the results presented hereafter will be based on the numerical solution of the Fokker–Planck equation (3.5) with $C_{I}\neq 0$ , where $\dot{\boldsymbol{p}}$ is given by the Jeffery equation for slender bodies and a description of the numerical method is detailed in Férec et al. (Reference Férec, Heniche, Heuzey, Ausias and Carreau2008). Hence, once the probability distribution function is numerically computed, the components of $\unicode[STIX]{x1D656}_{2}^{(m)}$ and $\unicode[STIX]{x1D656}_{4}^{(m)}$ are straightforwardly evaluated and the stress tensor components obtained.

3.2 Model predictions

In order to investigate the model predictions, a first set of simulations was performed imposing a fixed isotropic orientation distribution (i.e. ${\it\psi}=1/4{\rm\pi}$ ). All the results presented in this section are obtained for the following parameters: $K=50~\text{Pa}~\text{s}^{m}$ , ${\it\phi}=10\,\%$ , $a_{r}=20$ ( $L=280~{\rm\mu}\text{m}$ and $R=14~{\rm\mu}\text{m}$ ) and $h=D\sqrt{{\rm\pi}/4{\it\phi}}$ (Chung & Kwon Reference Chung and Kwon1996). This last parameter corresponds to the average distance between rods for an aligned orientation state (Dinh & Armstrong Reference Dinh and Armstrong1984) and is kept constant over the different simulations even if the particle microstructure orientation evolves.

A steady-state shear flow is considered (the imposed shear rate is $1~\text{s}^{-1}$ ), where subscripts 1, 2 and 3 stand for flow, velocity gradient and vorticity directions, respectively. Assuming an isotropic rod orientation distribution ( ${\it\psi}=1/4{\rm\pi}$ ), figure 2 reports the particle shear-stress contribution, ${\it\tau}_{12}^{p}$ , and the non-zero components of $\unicode[STIX]{x1D656}_{2}^{(m)}$ as functions of $m$ . ${\it\tau}_{12}^{p}$ exhibits a decrease with increasing shear-thinning character of the matrix (i.e. $m$ decreases). Starting from an isotropic effective rod orientation $\unicode[STIX]{x1D656}_{2}^{(m)}$ for $m=1$ , observation of $\unicode[STIX]{x1D656}_{2}^{(m)}$ reveals that the enhancement of the shear-thinning character of the matrix results in a system which behaves effectively like a suspension in which more and more rods would orient towards the vorticity axis: $\unicode[STIX]{x1D622}_{33}^{(m)}$ tends to 1 whereas $\unicode[STIX]{x1D622}_{22}^{(m)}$ and $\unicode[STIX]{x1D622}_{11}^{(m)}$ decrease to 0. Therefore, the effect of the scalar $|\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}$ is similar to a stretching force on the orientation distribution in the 3-direction as $\unicode[STIX]{x1D622}_{11}^{(m)}$ remains equal to $\unicode[STIX]{x1D622}_{22}^{(m)}$ .

Figure 2. Effect of the power-law index, $m$ , on the shear stress ${\it\tau}_{12}^{p}$ and on the non-zero components of $\unicode[STIX]{x1D656}_{2}^{(m)}$ .

The transient results are not shown here as the orientation states predicted by (3.5) are similar to the ones considering a Newtonian matrix and are well known (Férec et al. Reference Férec, Heniche, Heuzey, Ausias and Carreau2008). Note also that the influence of the rod aspect ratio and the particle concentration on the coupling coefficient, ${\it\mu}_{1}^{(m)}$ , have already been tested in Souloumiac & Vincent (Reference Souloumiac and Vincent1998).

3.3 Newtonian versus power-law fluids

Regarding the stress expressions for the rod-filled systems in both Newtonian (2.10) and power-law (3.2) matrices, some questions dealing with the enhancement/decrease of the suspension shear viscosity in a power-law matrix and the departure from Newtonian behaviour arise. In order to provide some answers, the specific shear viscosities of the suspensions ( ${\it\eta}_{sp}=({\it\eta}/{\it\eta}^{m})-1$ , where ${\it\eta}$ and ${\it\eta}^{m}$ are the shear viscosities for the composite and the unfilled matrix, respectively) in both Newtonian and pseudo-plastic matrices are given respectively by

(3.7) $$\begin{eqnarray}{\it\eta}_{sp}^{N}={\it\phi}\frac{2a_{r}^{2}}{3\log (h/R)}\langle p_{1}p_{1}p_{2}p_{2}\rangle ,\end{eqnarray}$$

and

(3.8) $$\begin{eqnarray}{\it\eta}_{sp}^{PL}={\it\phi}\frac{2a_{r}^{m+1}}{2^{m-1}(m+2)}\left[\frac{1-m}{m[1-(R/h)^{(1-m)/m}]}\right]^{m}\langle |\tilde{\dot{{\bf\gamma}}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}p_{1}p_{1}p_{2}p_{2}\rangle ,\end{eqnarray}$$

where $\tilde{\dot{{\bf\gamma}}}$ is the dimensionless strain-rate tensor. The quantities $\langle p_{1}p_{1}p_{2}p_{2}\rangle$ and $\langle |\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}p_{1}p_{1}p_{2}p_{2}\rangle$ are calculated from steady-state solutions of (3.5).

Figure 3 reports the Newtonian specific viscosity, ${\it\eta}_{sp}^{N}$ , as a function of Péclet number $(P_{e}=1/C_{I})$ , as well as ${\it\eta}_{sp}^{PL}$ for $m=0.6$ and 0.3, respectively (note that $C_{I}$ in (3.5) controls the diffusion due to rod interactions). Starting with an isotropic rod orientation state ( $P_{e}\ll 1$ ), the specific viscosity decreases with diminishing the power-law index, the largest value being observed for the suspension in a Newtonian fluid. In the region of $P_{e}\approx 10$ , the isotropic orientation distribution of the rods is destroyed and the resulting microstructure leads to systems having the most flow resistance. For Newtonian fluids, the maximum is reached for a configuration in which the largest number of rods is in close alignment with the principal axis of strain. At larger $P_{e}$ , decrease in viscosity due to rod orientation is observed and is more pronounced as $m$ becomes close to 1: rods orient toward directions in which the energy dissipation is the lowest. Finally at $P_{e}>10^{3}$ , the specific viscosity for the most shear-thinning matrix (i.e. $m=0.3$ ) becomes the highest.

Figure 3. Specific viscosities for rod suspensions in a Newtonian fluid ( $N$ ) and two power-law fluids with $m=0.6$ and 0.3, respectively.

In order to explain this behaviour, the $\unicode[STIX]{x1D622}_{11}^{(m)}$ and $\unicode[STIX]{x1D622}_{33}^{(m)}$ components of the second-order orientation tensor are depicted as functions of $P_{e}$ for the three different fluids in figure 4. As previously mentioned, considering an isotropic orientation state (i.e.  ${\it\psi}=\text{const.}$ ) for rods suspended in shear-thinning fluids (i.e.  $m\not =1$ ) leads to $\unicode[STIX]{x1D622}_{11}^{(m)}\not =\unicode[STIX]{x1D622}_{33}^{(m)}\not =1/3$ . At large $P_{e}$ and independently of $m$ , $\unicode[STIX]{x1D622}_{11}^{(m)}$ and $\unicode[STIX]{x1D622}_{33}^{(m)}$ tend to reach the same steady-state values close to 0.9 and 0.1, respectively. These results suggest that the shear-thinning behaviour of the matrix has no effect on the rod orientation dynamics in the case of strong flows. For $P_{e}\rightarrow \infty$ , the steady-state solution for the distribution function is simply a Dirac function that is ${\it\psi}(\boldsymbol{p})={\it\delta}(\boldsymbol{p}-{\bf\delta}_{1})$ , where the subscript 1 stands for the flow direction. Under these conditions, the steady-state shear viscosity for the suspension is equal to the unfilled viscosity because of the alignment of the rods in the planes of shear and the lack of inclusion of rod thickness in the model. Hence, as $P_{e}\rightarrow \infty$ , the specific viscosities converge to zero and all the components of $\unicode[STIX]{x1D656}_{2}^{(m)}$ are null, except $\unicode[STIX]{x1D622}_{11}^{(m)}$ , which is 1.

Figure 4. $\unicode[STIX]{x1D622}_{11}^{(m)}$ and $\unicode[STIX]{x1D622}_{33}^{(m)}$ as a function of $P_{e}$ for a Newtonian ( $N$ ) and two shear-thinning ( $m=0.6$ and 0.3) fluids, respectively.

4.1 Bi-viscous model

To investigate the departure from the plateau region, the rod contribution to the stress in the suspension is assumed to follow the Newtonian behaviour predicted by (2.15) at low shear rates and the power-law form given by (3.3) at high shear rates. The transition behaviour at intermediate shear rates is not necessary. Therefore, equating the shear viscosities of the suspension in both these limiting regions enables defining a critical shear rate $\dot{{\it\gamma}}_{c}^{{\it\phi}}$ at which the transition occurs from the zero-shear-rate plateau to the power-law regime (figure 5). After some straightforward calculus, the critical shear rate of the suspension is given by

(4.1) $$\begin{eqnarray}\dot{{\it\gamma}}_{c}^{{\it\phi}}=\left[\frac{{\it\eta}_{0}+{\it\eta}_{0}{\it\phi}\displaystyle \frac{2a_{r}^{2}}{3\log (h/R)}\langle p_{1}p_{1}p_{2}p_{2}\rangle }{K+K{\it\phi}\displaystyle \frac{2a_{r}^{m+1}}{2^{m-1}(m+2)}\left(\displaystyle \frac{1-m}{m[1-(R/h)^{(1-m)/m}]}\right)^{m}\langle |\tilde{\dot{{\bf\gamma}}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}|^{m-1}p_{1}p_{1}p_{2}p_{2}\rangle }\right]^{1/(m-1)}.\end{eqnarray}$$

The critical shear rate for an unfilled system, $\dot{{\it\gamma}}_{c}^{{\it\phi}=0}$ , is obtained by imposing ${\it\phi}=0$ in (4.1) (see also figure 5). The ratio ${\it\chi}=\dot{{\it\gamma}}_{c}^{{\it\phi}}/\dot{{\it\gamma}}_{c}^{{\it\phi}=0}$ (or inversely the ratio between the characteristic time for the unfilled and filled systems) is then introduced to investigate the departure from Newtonian behaviour. Hence, if ${\it\chi}<1$ , the rod suspension shear thins before the unfilled matrix and the opposite behaviour is observed for ${\it\chi}>1$ . Figure 6 depicts the ratio ${\it\chi}$ as a function of the Péclet number for three different values of the power-law index ( $m=0.9$ , 0.6 and 0.3, respectively). These results are obtained with ${\it\eta}_{0}=1000~\text{Pa}~\text{s}$ , ${\it\phi}=10\,\%$ and $K=1000~\text{Pa}~\text{s}^{m}$ and the dotted line in figure 6 indicates cases when the onset of rod suspension and unfilled matrix shear thinning occur at the same characteristic time. In the following analysis it will be useful to keep in mind the results for rod orientations presented in figure 4. The rod suspension shear thins before the unfilled matrix up to $\unicode[STIX]{x1D622}_{11}^{(m)}\leqslant 0.7$ ( $P_{e}\approx 10^{2}$ ), otherwise the opposite behaviour is observed and the range within which this phenomenon is noticed increases when decreasing the power-law index. Moreover, the largest departures from ${\it\chi}=1$ appear to match the specific viscosity peaks (see figure 3), where the average rod orientation is the closest to the principal axis of strain ( $P_{e}\approx 10$ ). As ${\it\eta}_{0}$ , $K$ and ${\it\phi}$ have been fixed arbitrary, curves plotted in figure 6 cannot be considered as master curves.

Figure 6. ${\it\chi}$ as a function of $P_{e}$ for three different values of the power-law index; $m=0.9$ , 0.6 and 0.3, respectively ( ${\it\phi}=10\,\%$ ). The dotted line indicates cases when the onset of rod suspension and unfilled matrix shear thinning occur at the same characteristic time.

4.2 Ellis model

The Ellis model has the advantage of expressing the shear rate as a function of the shear stress. In simple shear flow, the Ellis fluid response can then be written in the following form

(4.2) $$\begin{eqnarray}\frac{\dot{{\it\gamma}}}{{\it\tau}_{rz}}=\frac{1+(|{\it\tau}_{rz}|/{\it\tau}_{1/2})^{{\it\alpha}-1}}{{\it\eta}_{0}}.\end{eqnarray}$$

This model exhibits a Newtonian fluid character at ${\it\tau}_{rz}\rightarrow 0$ ( ${\it\eta}_{0}$ plateau viscosity) and a non-Newtonian power-law fluid character at high shear rates ( ${\it\alpha}$ is related to the power-law index by ${\it\alpha}=1/m$ ). While the index ${\it\alpha}$ is a measure of the degree of shear-thinning behaviour (the greater the value of ${\it\alpha}$ , the greater is the extent of shear thinning), ${\it\tau}_{1/2}$ represents the shear stress value at which the apparent matrix viscosity has dropped to half its zero-shear-rate value (at high shear stress or shear rate, $K=({\it\eta}_{0}{\it\tau}_{1/2}^{{\it\alpha}-1})^{1/{\it\alpha}}$ ). Therefore, the choice of an Ellis fluid appears to be a reasonable compromise as a constitutive model which captures more details of the viscosity trend observed experimentally. Substituting (4.2) into (2.7) leads to the following expression for the drag force

(4.3) $$\begin{eqnarray}\frac{z}{2}\dot{{\bf\gamma}}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}=\frac{f}{2{\rm\pi}{\it\eta}_{0}}\log (h/R)+\frac{f}{2{\rm\pi}{\it\eta}_{0}}\frac{1}{{\it\alpha}-1}\left(\frac{|f|}{2{\rm\pi}R{\it\tau}_{1/2}}\right)^{{\it\alpha}-1}[1-(R/h)^{{\it\alpha}-1}].\end{eqnarray}$$

The drag force $f$ has to be obtained numerically from (4.3) by means of simple numerical solvers based for example on a bi-section method. These do not require high computational time and generally do not restrict the applicability of the model. Hence, after computing the drag force numerically, the rod contribution to the stress tensor components are straightforwardly obtained using (2.8).

The model predictions are analysed by comparing them with experimental data obtained for glass-fibre-reinforced polystyrene melts. The experimental results are taken from the work of Chan et al. (Reference Chan, White and Oyanagi1978) and polystyrene (PS) polymers are filled respectively with 0 wt% (PS0), 20 wt% (PS20) and 40 wt% (PS40) of glass fibre of diameter equal to $12.7~{\rm\mu}\text{m}$ and length of order 2 mm for the PS20 and 1 mm for the PS40. A nonlinear least-squares algorithm is used to determine the Ellis parameters for the pure PS matrix and the parameters ${\it\eta}_{0}=750~\text{kPa}~\text{s}$ , ${\it\tau}_{1/2}=120~\text{kPa}$ and ${\it\alpha}=2.93$ were found. It can be seen in figure 7 that the viscosity of the unfilled polystyrene (PS0) is well predicted for shear rates spanning over six orders of magnitude.

Assuming a polystyrene density of ${\it\rho}_{PS}=1.00~\text{g}~\text{cm}^{-3}$ and a glass-fibre density of ${\it\rho}_{fib}=2.50~\text{g}~\text{cm}^{-3}$ , the calculated volume concentrations are found to be ${\it\phi}_{PS20}=9.1\,\%$ and ${\it\phi}_{PS40}=21.1\,\%$ for PS20 and PS40, respectively. The same fitting procedure is applied to obtain the shear viscosities of PS20 and PS40, as shown in figure 7. Note that the rod orientation state is given by the steady-state solution of the Fokker–Planck equation (3.5) (that is a non-isotropic orientation for the distribution function) and, $h$ appearing in (4.3) with $P_{e}$ are chosen as the two fitting parameters. Due to the important particle length, fibre breakage during the measurements cannot be ruled out and the values of $h$ and $P_{e}$ remain uncertain. The nonlinear least-squares algorithm gives $h=14.0~{\rm\mu}\text{m}$ and $P_{e}=8.70\times 10^{6}$ for PS20 and $h=7.8~{\rm\mu}\text{m}$ and $P_{e}=9.85\times 10^{5}$ for PS40. As can be seen in figure 7, the model predictions are in good agreement with the experimental findings. The obtained fitting values for $P_{e}$ suggest that the rod orientation state is in close to perfect alignment leading to a low contribution to the shear stress. However, as the coupling coefficients are due to the high aspect ratios (i.e. ${\approx}$ 160 and 80 for PS20 and PS40, respectively), the short glass fibres contribution to the shear stress cannot be omitted. We speculate that the measurement for the fibre dimensions have been performed before preparing the samples for the rheological tests, leading to lower aspect ratios due to fibre length breakage. Chan et al. (Reference Chan, White and Oyanagi1978) also reported the existence of normal stresses for the polystyrene systems but their measurements are not employed in our fitting procedure since the Ellis model predicts no normal stresses in shear flows for the unfilled matrix.

${\it\chi}$ can be estimated graphically by noting that $\dot{{\it\gamma}}_{c}^{{\it\phi}}=({\it\eta}_{PS20}^{N}|_{\dot{{\it\gamma}}=0.001\text{s}^{-1}}/{\it\eta}_{PS20}^{PL}|_{\dot{{\it\gamma}}=1\text{s}^{-1}})^{1/(m-1)}$ and $\dot{{\it\gamma}}_{c}^{{\it\phi}=0}=({\it\eta}_{PS0}^{N}|_{\dot{{\it\gamma}}=0.001\text{s}^{-1}}/{\it\eta}_{PS0}^{PL}|_{\dot{{\it\gamma}}=1\text{s}^{-1}})^{1/(m-1)}$ , where ${\it\eta}_{PS0}^{N}|_{\dot{{\it\gamma}}=0.001\text{s}^{-1}}\approx 7.5\times 10^{5}~\text{Pa}~\text{s}$ , ${\it\eta}_{PS20}^{PL}|_{\dot{{\it\gamma}}=1\text{s}^{-1}}\approx 5.0\times 10^{5}~\text{Pa}~\text{s}$ , ${\it\eta}_{PS0}^{N}|_{\dot{{\it\gamma}}=0.001\text{s}^{-1}}\approx 7.5\times 10^{5}~\text{Pa}~\text{s}$ and ${\it\eta}_{PS0}^{PL}|_{\dot{{\it\gamma}}=1\text{s}^{-1}}\approx 2.5\times 10^{5}~\text{Pa}~\text{s}$ , respectively. Note that ${\it\eta}_{PS0}^{PL}|_{\dot{{\it\gamma}}=1\text{s}^{-1}}$ and ${\it\eta}_{PS20}^{PL}|_{\dot{{\it\gamma}}=1\text{s}^{-1}}$ may be obtained by extrapolation of the power-law viscosity curves up to $\dot{{\it\gamma}}=1\text{s}^{-1}$ (in our case, this extrapolation is not necessary). It is then found that ${\it\chi}\approx 0.6<1$ , indicating that the rod suspension shear thins before the unfilled matrix, as also observed experimentally in figure 7.

Figure 7. Comparison of the steady-shear viscosity for the glass-fibre-filled polystyrene systems (Chan et al. Reference Chan, White and Oyanagi1978) with the model predictions.

4.3 Carreau model

Based on the previous results, we can extend our analysis to rod suspensions in Carreau fluids. By neglecting the viscosity plateau at large shear rates, the Carreau model writes

(4.4) $$\begin{eqnarray}{\it\eta}={\it\eta}_{0}[1+({\it\lambda}_{0}\dot{{\it\gamma}})^{2}]^{(m-1)/2},\end{eqnarray}$$

where the parameter ${\it\lambda}_{0}$ is a time constant for the fluid. The value of ${\it\lambda}_{0}$ determines the shear rate at which the transition occurs from the zero-shear-rate plateau to the power-law portion. This time, the micromechanical analysis at the rod scale is numerically performed: equation (2.5) is solved along incremental particle segments to express the shear stress at the wall considering that the fluid follows the Carreau model. Hence, the extra stress due to the particle contribution (2.8) is numerically computed and the material functions for the suspensions can be deduced.

Figure 8 depicts the specific shear viscosities, ${\it\eta}_{sp}$ , when the rods are suspended in a Newtonian, power-law and Carreau fluids, respectively. For the following simulation results, an isotropic rod orientation distribution is assumed ( ${\it\psi}=1/4{\rm\pi}$ ), the Newtonian viscosity is ${\it\eta}_{0}=1000~\text{Pa}~\text{s}$ , the consistency is given by $K={\it\eta}_{0}{\it\lambda}_{0}^{m-1}$ , where the power-law index is $m=0.3$ and the characteristic time for the Carreau model is ${\it\lambda}_{0}=0.1~\text{s}$ . The particle phase properties are ${\it\phi}=10\,\%$ , $a_{r}=20$ and $h=D\sqrt{{\rm\pi}/4{\it\phi}}$ . As expected, at low shear rates, the specific viscosity obtained for rods suspended in a Carreau fluid matched the one predicted for a Newtonian suspending fluid (2.15) and at high shear rates, the observed trend is the power-law behaviour (3.3).

This approach to predicting the material functions for rod suspensions in Carreau fluid is useful for viscometric flows but leads to the high computational times required for non-homogenous problems. Although a micromechanical description at the particle scale is of key importance to the understanding of the flow kinematics of suspensions in shear-thinning fluids, the approach used in the case of a Carreau fluid appears to be rather restricted for practical applications, at least for the moment.

Figure 8. Specific viscosities for rod suspensions in Newtonian, power-law and Carreau fluids, respectively.