# Rheological Modeling and Forming Process Simulation of CNT Nanocomposites

Elías Cueto, Rosa Monge, Francisco Chinesta, Arnaud Poitou, Iciar Alfaro, Malcolm Mackley

► To cite this version:

Elías Cueto, Rosa Monge, Francisco Chinesta, Arnaud Poitou, Iciar Alfaro, et al.. Rheological Mod-
eling and Forming Process Simulation of CNT Nanocomposites. International Journal of Material
Forming, 2010, 3 (2), pp.1327-1338. 10.1007/s12289-009-0659-6 . hal-01007260

HAL Id: hal-01007260

https://hal.science/hal-01007260

Submitted on 3 Feb 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Rheological modeling and forming process simulation of CNT nanocomposites

Elias Cueto · Rosa Monge · Francisco Chinesta · Arnaud Poitou · Iciar Alfaro · Malcolm R. Mackley

E. Cueto · R. Monge · I. Alfaro
I3A, Universidad de Zaragoza, Maria de Luna
7, Zaragoza 50012, Spain
e mail: ecueto@unizar.es

F. Chinesta · A. Poitou
GEM UMR CNRS Centrale Nantes,
1 rue de la Noe, BP 92101, 44321 Nantes cedex 3, France

M. R. Mackley
Department of Chemical Engineering, University of Cambridge,
Pembroke Street,
Cambridge CB2 3RA, UK

**Abstract** This paper is concerned with the rheological modelling of untreated Carbon Nanotubes (CNTs) suspended within an epoxy resin and its application in nanocomposites forming processes. The untreated CNT suspensions exhibited significant steady shear-thinning and contained optically resolvable aggregate structures. A simple orientation model, based on a Fokker-Planck advection-diffusion description, failed to capture the experimentally observed rheological responses for untreated CNT suspensions. A new model named the “Aggregation/Orientations” (AO) model has been developed to describe the experimental findings. The model integrated elements of both a standard orientation model and aggregation modelling concepts within the Fokker-Planck formalism. However, the proposed model has different weakness; the most important concerns the affine orientation of tubes with the fluid deformation independently of the state of aggregation. This mechanism seems to be in opposition to the experimental observations. In this paper we propose an alternative rheological model that is then employed for simulating film processing by employing spin coating.

**Keywords** Carbon nanotube suspensions · Kinetic theory · Meshless methods · Natural element method · Spin coating

## Introduction

The experimental rheology of Carbon Nanotubes (CNT) has been studied extensively over the last five years (see for example, [11, 18, 20]) and the rheological behaviour of treated and untreated CNT suspensions has been found to be significantly different particularly in relation to shear thinning.

In terms of direct visualisation, there are publications on the aggregates or bundles of CNTs [9, 12, 13, 20]. In general, optical microstructures were observed in untreated CNTs, whereas suspensions with treated CNTs tended to be more optically transparent [15, 16, 22], making structural characterisation of treated CNT suspensions difficult.

In terms of modelling, the evolution of flow curves for CNT suspensions can be described by different empirical models such as the Cross model, the Carreau model and the Krieger-Dougherty model. A stress-microstructure model, however, is desirable as information such as CNT orientation and state of aggregation is intimately related to the final physical properties of CNT composite materials. The current CNT modelling is based on both fibre orientation models and transient network modelling concepts.

In some of our former works we proposed an orientation based model for describing the rheology of functionalized CNT [17]. That orientation based model was applied for simulating spin coating processes involving semi-dilute suspensions of functionalized CNTs [4]. However, this model was not able to describe the rheological behaviour of non-functionalized CNT suspensions. In that case an orientation/aggregation based model was proposed to take into account the presence and evolution of aggregates whose impact on the suspension rheology is significant (see Ma et al. [14]). However, the orientation/aggregation model proposed in Ma et al. [14] has some weakness because the
tube orientation was assumed governed by the Jeffery equation that assuming the aspect ratio of tubes predicts a quasi-affine orientation of tubes with the matrix fluid deformation. This prediction seems to be in contradiction with some experimental observations that indicate that tubes belonging to large aggregates do not orient as the Jeffery's model predicts.

The first part of the paper revisits the orientation/aggregation model proposed in Ma et al. [14]. Then, the question related to the delayed orientation is deeply addressed. A two populations model similar to the one considered for modelling associative polymers [8, 21] is then developed. In that model the orientation velocity is different for both populations. Finally the model is applied for simulating spin coating.

## Experiments

When dispersed in an epoxy resin, both treated and untreated CNTs showed steady shear-thinning character-

istics, but to a different extent (Fig. 1). All the details concerning the experimental protocol can be found in Ma et al. [14, 17, 16]. The epoxy matrix is optically transparent, essentially Newtonian and with a base viscosity of 10 Pa.s. As shown in the figure, untreated CNT suspensions possessed a higher apparent low-shear viscosity coupled with optical microstructure of CNT aggregates, whereas treated CNT suspensions showed no optically resolvable CNT aggregates with a less pronounced viscosity enhancement effect.

The steady shear rheology of untreated CNT suspensions was studied in greater detail and the evolution of apparent shear viscosity ($\eta_a$) as a function of shear rate is shown in Fig. 2 for four different concentrations. The untreated CNT suspensions exhibited a significant viscosity enhancement at low shear rates and it was observed that the addition of only 0.1% CNT had led to an order-of-magnitude increase in $\eta_a$. At high shear rates, the viscosity of the suspensions decreased asymptotically to that of the base epoxy and these findings are consistent with those reported by Rahatekar et al. [20].

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Apparent shear viscosity ($\eta_a$) of **a** epoxy, **b** 0.25% treated CNT and **c** 0.25% untreated CNT suspensions. Treated and untreated CNTs both exhibited shear thinning characteristics but to a different extent when dispersed in an epoxy matrix. Untreated CNTs showed a higher low shear viscosity coupled with optical microstructures of CNT</figcaption></figure>

aggregates whereas treated CNTs showed no clear optically resolvable CNT aggregates with a less pronounced viscosity enhancement effect. (Treated CNTs were supplied by Nanocomposites Inc, USA and the same mixing protocol was used for preparing both treated and untreated CNT suspensions.)
<figure><img src="image_3.png" /><figcaption>Fig. 2 Apparent shear viscosity (η_a) as a function of shear rate for a 0.25% CNT, b 0.1%, c 0.05%, d 0.025% CNT suspensions and e epoxy only. The η_a data presented are experimental data and the values were obtained after a shearing time of 100 s</figcaption></figure>

In terms of microstructure, untreated CNT suspensions showed a clear optical texture which evolved depending on the shear conditions (see for example Fig. 3a–d). At low shear rates, coarsening of the optical texture was generally observed after applying steady shear for a certain period of time (as shown in Fig. 3b). At high shear rates, CNT suspensions were more uniformly dispersed and the aggregate size of CNTs decreased as shear rate increased (Fig. 3c and d).

## Rheological modelling

### The orientation model

In the context of short-fibre suspensions modelling, an orientation model can be used to describe the evolution of steady shear viscosity by coupling the flow kinematics with fibre orientation. First, an orientation model was assumed where the CNTs were considered as short, rigid fibres that can rotate and align in a shear flow and where the evolution of the viscosity contribution due to the presence of CNTs depended only on the orientation of CNTs.

The momentum balance equation (neglecting the inertia and mass terms):

$$ \nabla \cdot \sigma = \mathbf{0} \qquad (1) $$

where $\sigma$ is the stress tensor. The mass balance equation for an incompressible fluid reads:

$$ \nabla \cdot \mathbf{v} = \mathbf{0} \qquad (2) $$

where $\mathbf{v}$ represents the velocity field.

The constitutive equation for a dilute suspension of high aspect-ratio particles (with negligible Brownian diffusion) is usually written as:

$$ \sigma = -P\mathbf{I} + 2\eta(\mathbf{D} + N_p(\mathbf{a}_4 : \mathbf{D})) = -P\mathbf{I} + 2\eta\mathbf{D} + \tau \quad (3) $$

where $P$ denotes the hydrostatic pressure, $\mathbf{I}$ is the unit tensor, $\eta$ is the matrix viscosity, $\mathbf{D}$ is the strain-rate tensor, $N_p$ is a scalar parameter that depends on fiber concentration and aspect ratio, “:” is the tensor product twice contracted (i.e. $(\mathbf{a}_4 : \mathbf{D})_{ij} = (\mathbf{a}_4)_{ijkl}\mathbf{D}_{kl}$) and $\mathbf{a}_4$ is the fourth order orientation tensor and it is defined as:

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

where $\mathbf{p}$ is the unit vector aligned in the fiber axis direction, “⊗” denotes the tensorial product (i.e. $(\mathbf{p} \otimes \mathbf{p})_{ij} = \mathbf{p}_i \mathbf{p}_j$), and $\psi(\mathbf{p})$ is the orientation distribution function which satisfies the normality condition $\int \psi(\mathbf{p}) d\mathbf{p} = 1$. In general, the distribution function depends on the physical coordinates $\mathbf{x}$, the time $t$ and the conformation coordinates $\mathbf{p}$ (i.e. $\psi(\mathbf{x},t,\mathbf{p})$), but if the flow is homogeneous and steady, the fibre orientation distribution would only depend on the conformation coordinates, (i.e. $\psi(\mathbf{p})$). For dilute and non-aggregating short fibre suspensions, different analytical expressions for the parameter $N_p$ are available [19]. However, those expressions do not immediately apply to semi-dilute/concentrated or aggregating suspensions. Aggregation has been observed within the semi-dilute untreated CNT suspensions studied in this paper and $N_p$ is therefore determined from experimental data fitting.

If all the fibres are completely aligned in the direction the fourth order orientation tensor can be written in the following form:

$$ \mathbf{a}_4 = \mathbf{a}_2 \otimes \mathbf{a}_2 \qquad (5) $$

where $\mathbf{a}_2 = \int \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}) d\mathbf{p}$ denotes the second order orientation tensor and corresponding components are given as $(\mathbf{a}_4)_{ijkl} = (\mathbf{a}_2)_{ij}(\mathbf{a}_2)_{kl}$.

In cases where the fibres are not completely aligned in one direction, Eq. (5) is not exact, but it provides an approximation for the fourth-order orientation tensor. Equation (5) is also known as the quadratic closure relation: $\mathbf{a}_4^{quad} = \mathbf{a}_2 \otimes \mathbf{a}_2$ and other forms of closure relations are also available in the literature (see for example Advani and Tucker [1]).

The orientation distribution function $\psi(\mathbf{p})$ can be determined by solving its balance equation, which is also known as the Fokker-Planck equation:

$$ \frac{d\psi(\mathbf{p})}{dt} + \frac{\partial}{\partial \mathbf{p}} \left( \psi(\mathbf{p}) \frac{d\mathbf{p}}{dt} \right) = 0 \qquad (6) $$

where $\frac{d\psi}{dt}$ is the material derivative: $\frac{d\psi}{dt} = \frac{\partial\psi}{\partial t} + (\mathbf{v} \cdot \nabla)\psi$.
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Optical microstructure of 0.1% CNT suspended within an epoxy resin after shearing at **a** 0 s⁻¹, **b** 0.5 s⁻¹, **c** 10 s⁻¹ and **d** 100 s⁻¹. The micrographs were captured after a shearing time of 100 s and using the Cambridge Shear System (CSS 450) with an optical depth of 130 µm</figcaption></figure>

Once the distribution function is calculated, the computation of its moments (orientation tensors) then becomes straightforward and does not require any closure approximation.

For a dilute suspension of spheroids and in the absence of Brownian motion, the evolution of fibre orientation was described by Jeffery [10]:

$$ \frac{d\mathbf{p}}{dt} = \mathbf{\Omega} \cdot \mathbf{p} + \lambda (\mathbf{D} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p})\mathbf{p}) \quad (7) $$

where $\mathbf{\Omega}$ is the vorticity tensor and $\lambda$ is a constant that depends on the fibre aspect ratio $r$ ($r$ = fiber length / fibre diameter) and $\lambda$ is given as:

$$ \lambda = (r^2 - 1)/(r^2 + 1) \quad (8) $$

For particles with an infinite aspect ratio $\lambda \approx 1$, Eq. (7) reduces to:

$$ \frac{d\mathbf{p}}{dt} = \nabla\mathbf{v} \cdot \mathbf{p} - (\mathbf{p}^T \cdot \mathbf{D} \cdot \mathbf{p})\mathbf{p} \quad (9) $$

By introducing Eq. (7) into the expression defining the second order orientation tensor, the equation governing the evolution of the second order orientation tensor can be deduced:

$$ \frac{d\mathbf{a}_2}{dt} = \mathbf{\Omega} \cdot \mathbf{a}_2 - \mathbf{a}_2 \cdot \mathbf{\Omega} + \lambda (\mathbf{D} \cdot \mathbf{a}_2 + \mathbf{a}_2 \cdot \mathbf{D} - 2(\mathbf{a}_4 : \mathbf{D})) \quad (10) $$

In the case of semi-dilute fibre suspensions, Folgar and Tucker [6] proposed the introduction of a diffusion term

$(D_r)$ in Eq. (6) to account for fiber interaction effects, resulting in the following expression:

$$ \frac{d\psi(\mathbf{p})}{dt} + \frac{\partial}{\partial \mathbf{p}} \{\psi(\mathbf{p}) \frac{d\mathbf{p}}{dt}\} = \frac{\partial}{\partial \mathbf{p}} \{D_r \frac{\partial \psi(\mathbf{p})}{\partial \mathbf{p}}\} \quad (11) $$

The diffusion coefficient $D_r$ takes into account the Brownian effects and all the randomizing effects that could be originated by the tubes interactions. In the present work this rotary diffusion is assumed constant but its dependence with the shear rate could be introduced without major difficulties.

The equation governing the evolution of $\mathbf{a}_2$ (in 3-D) then becomes:

$$ \frac{d\mathbf{a}_2}{dt} = \mathbf{\Omega} \cdot \mathbf{a}_2 - \mathbf{a}_2 \cdot \mathbf{\Omega} + k(\mathbf{D} \cdot \mathbf{a}_2 + \mathbf{a}_2 \cdot \mathbf{D} - 2(\mathbf{a}_4 : \mathbf{D})) - 6D_r \left(\mathbf{a}_2 - \frac{\mathbf{I}}{3}\right) \quad (12) $$

Equation (12) shows that the evolution equation for the second-order orientation tensor involves the fourth-order orientation tensor and in order to solve the problem in a closed form, closure approximations which express the fourth-order orientation tensor in terms of lower-order tensors are commonly used [1, 5]. Instead of using closure approximations, a number of authors solved the Fokker-Planck description directly (see for example, [2, 3]). However, the latter approach also requires the computation of a multi-dimensional distribution function and thus involves considerable computation efforts.
In the Fokker-Planck based orientation model, $N_p$ and $D_r$ are the key adjustable parameters to fit with experimental data. The model was found to be successful in capturing the shear thinning characteristic for some surface-treated CNT suspensions [17], where no optically resolvable CNT aggregate was observed. The success of this model for treated CNTs can be understood in terms of the non-aggregating characteristic of treated CNT suspensions, where CNT orientation was expected to be the key factor in controlling the evolution of steady shear viscosity.

The Fokker-Planck based orientation model with the use of a single $D_r$ and $N_p$ values, however, failed to fit the experimental flow curves for untreated CNT suspensions, where the viscosity enhancement was very pronounced at low shear rates. It is believed that the extra contribution to shear viscosity at low shear rates is due to the presence of CNT aggregates, which were experimentally observed in optical microstructure studies (see also Fig. 3b-d).

### The orientation/aggregation model

The FP orientation model was unable to describe the shear-thinning characteristic for untreated CNT suspensions, and therefore, a new model named the Aggregation/Orientational

(AO) model was developed. The new model considered a hierarchy of CNT aggregate structures within an untreated CNT suspension, where the shear viscosity was controlled by the state of aggregation and CNT orientation. The Fokker-Planck description was modified to incorporate aggregation/disaggregation kinetics and the derivation of the AO model is summarized in this section.

The aggregation/orientation distribution function in the AO model is written as $\psi(\mathbf{x},t,\mathbf{p},n)$, where $n \in [0,1]$ describes the state of aggregation ($n=0$ corresponds to CNTs that are free from entanglement and $n=1$ represents a CNT aggregate network). It is assumed that the flow is steady and homogeneous and that different populations ($n$) are present within the control volume considered in the balance equation.

For an arbitrary population $n$ (where $n \neq 0$ nor $1$), the population can increase as a result of the aggregation of smaller aggregates ($r < n$) or the disaggregation of larger aggregates ($r > n$). On the other hand, the population $n$ can decrease because of the disaggregation of population $n$ forming less entangled aggregates ($r < n$) or the aggregation of $n$ forming more entangled aggregates ($r > n$). If a constant aggregation velocity ($v_c$) and a constant disaggregation velocity of ($v_d$) are assumed, the following balance expression can be written for the population $n$:

$$ \frac{\partial \psi(\mathbf{p}, n)}{\partial t} \bigg|_{\text{Aggregation}} \, dn \, d\mathbf{p} \, dt = \left\{ v_c \int_{r=0}^{r=n} \psi(\mathbf{p}, r) dr + v_d \int_{r=n}^{r=1} \psi(\mathbf{p}, r) dr - v_d \psi(\mathbf{p}, n) \int_{r=0}^{r=n} dr - v_c \psi(\mathbf{p}, n) \int_{r=n}^{r=1} dr \right\} dn \, d\mathbf{p} \, dt \quad (13) $$

The balance equation considering orientation only is given as:

$$ \frac{\partial \psi(\mathbf{p}, n)}{\partial t} \bigg|_{\text{Orientation}} \, dn \, d\mathbf{p} \, dt = \left\{ -\frac{\partial}{\partial \mathbf{p}} \left\{ \psi(\mathbf{p}, n) \frac{d\mathbf{p}}{dt} \right\} + \frac{\partial}{\partial \mathbf{p}} \left\{ D_r(n) \frac{\partial \psi(\mathbf{p}, n)}{\partial \mathbf{p}} \right\} \right\} dn \, d\mathbf{p} \, dt \quad (14) $$

The Fokker-Planck equation which includes both CNT orientation and aggregation/disaggregation kinetics becomes:

$$ \frac{\partial \psi(\mathbf{p}, n)}{\partial t} \bigg|_{\text{Orientation}} + \frac{\partial \psi(\mathbf{p}, n)}{\partial t} \bigg|_{\text{Aggregation}} \left\{ \frac{\partial}{\partial \mathbf{p}} \left\{ \psi(\mathbf{p}, n) \frac{\partial \mathbf{p}}{\partial t} \right\} + \frac{\partial}{\partial \mathbf{p}} \left\{ D_r(n) \frac{\partial \psi(\mathbf{p}, n)}{\partial t} \right\} \right\} + \left\{ v_c \int_{r=0}^{r=n} \psi(\mathbf{p}, r) dr + v_d \int_{r=n}^{r=1} \psi(\mathbf{p}, r) dr \quad v_d \psi(\mathbf{p}, n) \int_{r=0}^{r=n} dr \quad v_c \psi(\mathbf{p}, n) \int_{r=n}^{r=1} dr \right\} \quad (15) $$

For a homogeneous flow at steady state,

$$ \left\{ -\frac{\partial}{\partial \mathbf{p}} \left\{ \psi(\mathbf{p}, n) \frac{d\mathbf{p}}{dt} \right\} + \frac{\partial}{\partial \mathbf{p}} \left\{ D_r(n) \frac{\partial \psi(\mathbf{p}, n)}{\partial \mathbf{p}} \right\} \right\} + \left\{ v_c \int_{r=0}^{r=n} \psi(\mathbf{p}, r) dr + v_d \int_{r=n}^{r=1} \psi(\mathbf{p}, r) dr - v_d \psi(\mathbf{p}, n) \int_{r=0}^{r=n} dr - v_c \psi(\mathbf{p}, n) \int_{r=n}^{r=1} dr \right\} = 0 \quad (16) $$

The solution of Eq. (16) should verify the normality condition $\int_{S(0,1)} \int_{0}^{1} \psi(\mathbf{p}, n) dnd\mathbf{p} = 1$, where $S(0,1)$ represents the unit sphere surface centred at the origin of the orientation conformation space.
The distribution function $\psi(\mathbf{p},n)$ can be obtained by solving Eq. (16) and the fourth order orientation tensor can be expressed as:

$$ \mathbf{a}_4(n) = \int \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \psi(\mathbf{p}, n) d\mathbf{p} \qquad (17) $$

The stress contribution due to the presence of CNTs for a multi-population system can be obtained by generalising Eq. (3)

$$ \tau = 2\eta \int_{0}^{1} N_{p}(n)(\mathbf{a}_{4}(n) : \mathbf{D}) dn \qquad (18) $$

The rheological model described above provides a first approximation for describing a suspension with aggregating CNTs and contains four rheological functions: the concentration and aspect-ratio parameter $N_p$, the rotary diffusion coefficient ($D_r$) and the aggregation and disaggregation velocities ($v_c$ and $v_d$). To simplify the analysis and minimize the degrees of freedom, the following assumptions are made, but other more complex models could be proposed:

a. $N_p$ varies linearly with the population variable $n$. If we assume that $N_p^{\min} \ll N_p^{\max}$ and that $N_p^{\min} \approx 0$, then $N_p(n) \approx N_p^{\max}n$ (see Eq. (19a)). (Note: $N_p^{\min} = N_p(n=0)$ and $N_p^{\max} = N_p(n=1)$). This approximation is valid when aggregate populations are present, because in that case the contribution of free tubes can be neglected with respect to the aggregates contribution to the stress. Obviously, when the shear rate becomes too high and aggregates are broken, one should introduce a non-zero $N_p^{\min}$. In process conditions the presence of aggregates is always present justifying this crude approximation.

b. In the modelling of associative polymers, a number of authors [8, 21] assumed that chains within an aggregate network can undergo affine deformation and therefore $D_r(n=1) \approx 0$. In the case of CNT suspensions, there is however experimental evidence showing that the opposite is true. CNT aggregates tended to rotate in simple shear flow and CNTs within aggregate were observed to be more isotropically orientated compared with disentangled CNTs [9, 13]. Despite all these evidences in Ma et al. [14] the rotary diffusion coefficient was assumed to increase linearly with the population $n$ as given in Eq. (19b) to enforce an isotropic orientation inside the aggregates.

c. $v_d$ increases linearly with the shear rate and it takes the form as given in Eq. (19c) where $\dot{\gamma}_{\max}$ is a characteristic shear rate above which the suspension viscosity coincides with the matrix viscosity. At high shear rates, $v_d$ remains constant and $v_d = v_d^{\max}$ (19d).

d. $v_c$ decreases linearly with the shear rate and it takes the form as given in Eq. (19e) and for a shear rate higher than $\dot{\gamma}_{\max}$, the aggregation velocity becomes zero (Eq. (19f)).

$$ \left\{ \begin{array}{l} N_p(n) = N_p^{\max}n + N_p^{\min}(1-n) \approx N_p^{\max}n \quad (a) \\ D_r(n) = D_r^{\max}n \quad (b) \\ v_d(\dot{\gamma}) = \begin{cases} v_d^{\max} \left(\frac{\dot{\gamma}}{\dot{\gamma}_{\max}}\right) & \text{if } \dot{\gamma} \le \dot{\gamma}_{\max} \\ v_d^{\max} & \text{if } \dot{\gamma} > \dot{\gamma}_{\max} \end{cases} \quad (c) \\ v_c(\dot{\gamma}) = \begin{cases} v_c^{\max} \left(\frac{\dot{\gamma}_{\max} - \dot{\gamma}}{\dot{\gamma}_{\max}}\right) & \text{if } \dot{\gamma} \le \dot{\gamma}_{\max} \\ 0 & \text{if } \dot{\gamma} > \dot{\gamma}_{\max} \end{cases} \quad (e) \\ v_d^{\max} \quad (f) \end{array} \right. \qquad (19) $$

This model was successfully applied for describing the non-linear rheological behavior of non-functionalized CNTs suspensions. The model only involves 5 rheological parameters easily identified by fitting the experimental flow curves. However, its use in forming processes simulation, where complex flows geometries are encountered, remains quite delicate because the significant computing time cost associated to the solution of the transient orientation/aggregation kinetic theory model.

On the other hand, and from a conceptual point of view, we must accept that the hypothesis (b) just introduced is not founded in any physical reasoning. Simply, as we know that the diffusion coefficient tends to induce an isotropic tube distribution and according to the observations [9, 13] this is the orientation state observed in aggregates, we decided to increase artificially the diffusion coefficient to be able to predict the observed state of the microstructure.

## Delaying orientation

In this section we are trying to propose a mechanism that could explain the lower observed order parameter related to aggregates.

For this purpose we start given a physical interpretation of the Jeffery equation. We consider the case of a rigid tube, with negligible mass and with a quasi infinite aspect ratio, immersed in a Newtonian fluid flow characterized by the velocity field that is assumed non perturbed by the presence of the tube. The dimension of the tube being micrometric, the shear rate can be assumed constant in the scale of the tube. The length of the tube is $2L$, the origin of the system of coordinates is assumed located on the tube centre of gravity and the tube orientation is defined by a unit vector $\mathbf{p}$ aligned with the tube axis. The tube is idealized by using two dumbbells located at the tube ends and on which the different forces apply, joined by rigid connector of length $2L$ (see the schema in Fig. 4).
<figure><img src="image_5.png" /><figcaption>Fig. 4 Schema of a massless ($m \approx 0$) tube suspended within a flow</figcaption></figure>

We are assuming that the tube is subjected to hydrodynamic drag forces $\mathbf{F}_H$, Brownian forces being ignored by the moment. The hydrodynamic forces write:

$$ \mathbf{F}_H = \xi(\mathbf{v} - L\dot{\mathbf{p}}) \qquad (20) $$

where $\mathbf{v}$ is the fluid velocity at the position occupied by the dumbbell (as just commented this velocity is assumed undisturbed by the presence of the tube) and $L\dot{\mathbf{p}}$ represents the velocity of the dumbbell. Thus, Eq. (20) can be rewritten as:

$$ \mathbf{F}_H = \xi L(\nabla \mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}) \qquad (21) $$

The equilibrium of forces in ensured because if $\mathbf{F}_H$ applies on a dumbbell, $-\mathbf{F}_H$ applies in the opposite one. Now, the inertia moment vanishing, the resulting torque must also vanish, i.e.

$$ \mathbf{M} = 2L\mathbf{p} \times \mathbf{F}_H = \mathbf{0} \qquad (22) $$

which implies that the force and the tube are collinear, i.e.:

$$ \mathbf{F}_H = \alpha \mathbf{p} \qquad (23) $$

or

$$ \xi L(\nabla \mathbf{v} \cdot \mathbf{p} - \dot{\mathbf{p}}) = \alpha \mathbf{p} \qquad (24) $$

Multiplying both terms by $\mathbf{p}$ and taking into account the tube inextensibility ($\mathbf{p} \cdot \mathbf{p} = 1$ that also implies $\dot{\mathbf{p}} \cdot \mathbf{p} = 0$), it results:

$$ \alpha = \xi L(\mathbf{p} \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \qquad (25) $$

Introducing (25) into (23) we obtain the expression of the hydrodynamic force:

$$ \mathbf{F}_H = \xi L(\mathbf{p} \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \mathbf{p} \qquad (26) $$

that as expected is aligned in the tube direction. Concerning the orientation kinematics, it could be derived by introducing Eq. (25) into (24) that results in:

$$ \dot{\mathbf{p}} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p} \cdot \nabla \mathbf{v} \cdot \mathbf{p})\mathbf{p} = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p})\mathbf{p} \qquad (27) $$

that corresponds to the Jeffery equation for ellipsoids with infinite aspect ratio. This equation implies that tubes,

assumed as particles with infinite aspect ratio, orient affinely with the matrix fluid deformation.

For introducing a delay in the tube orientation with respect to the fluid deformation we decided pragmatically to introduce an additional force that scales with the dumbbell velocity and tends to resist its movement. This force that we denote by $\mathbf{F}_R$, writes:

$$ \mathbf{F}_R = -\beta(L\dot{\mathbf{p}}) \qquad (28) $$

The equilibrium of forces is not affected by the introduction of this additional force; however the equilibrium of torques must be enforced:

$$ \mathbf{M} = 2L\mathbf{p} \times (\mathbf{F}_H + \mathbf{F}_R) = \mathbf{0} \qquad (29) $$

that implies:

$$ \xi L \left( \nabla \mathbf{v} \cdot \mathbf{p} - \left( 1 + \frac{\beta}{\xi} \right) \dot{\mathbf{p}} \right) = \alpha \mathbf{p} \qquad (30) $$

Multiplying again both terms by $\mathbf{p}$, and taking into account the tube inextensibility, it results again:

$$ \alpha = \xi L(\mathbf{p} \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \qquad (31) $$

from which, it results combining Eqs. (30) and (31):

$$ \left(1 + \frac{\beta}{\xi}\right) \dot{\mathbf{p}} = \left(\nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p})\mathbf{p}\right) \qquad (32) $$

or

$$ \dot{\mathbf{p}} = \kappa(\nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p})\mathbf{p}) \qquad (33) $$

where $\kappa = \frac{\xi}{\xi+\beta}$. It is easy to notice that Eq. (33) corresponds to the Jeffery expression if a sliding factor is introduced into the fluid velocity:

$$ \dot{\mathbf{p}} = (\nabla \mathbf{v}' \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D}' \cdot \mathbf{p})\mathbf{p}) \qquad (34) $$

with $\mathbf{v}' = \kappa \mathbf{v}$. However, it is well known that this model that was proposed in the past does not fulfil the objectivity principle.

For avoiding the objectivity crime, we consider an objective force $\mathbf{F}_R$:

$$ \mathbf{F}_R = -\beta L(\dot{\mathbf{p}} - \Omega \cdot \mathbf{p}) \qquad (35) $$

Now, we follow again the previous procedure for deriving the tube rotational velocity:

$$ \xi L \left( \nabla \mathbf{v} \cdot \mathbf{p} - \left( 1 + \frac{\beta}{\xi} \right) \dot{\mathbf{p}} + \frac{\beta}{\xi} \Omega \cdot \mathbf{p} \right) = \alpha \mathbf{p} \qquad (36) $$

Multiplying again both terms by $\mathbf{p}$, and taking into account the tube inextensibility and the anti-symmetry of the vorticity tensor, it results again:

$$ \alpha = \xi L(\mathbf{p} \cdot \nabla \mathbf{v} \cdot \mathbf{p}) \qquad (37) $$
from which, it results combining Eqs. (36) and (37):

$$ (1 + \frac{\beta}{\xi}) \dot{\mathbf{p}} = \left( \nabla \mathbf{v} \cdot \mathbf{p} + \frac{\beta}{\xi} \Omega \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p}) \mathbf{p} \right) \quad (38) $$

Now, introducing the relationship $\nabla \mathbf{v} = \mathbf{D} + \Omega$ into Eq. (38) it results:

$$ (1 + \frac{\beta}{\xi}) \dot{\mathbf{p}} = \left( \left( 1 + \frac{\beta}{\xi} \right) \Omega \cdot \mathbf{p} + \mathbf{D} \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p}) \mathbf{p} \right) \quad (39) $$

or

$$ \dot{\mathbf{p}} = (\Omega \cdot \mathbf{p} + \kappa(\mathbf{D} \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p})\mathbf{p})) \quad (40) $$

that corresponds to the usual Jeffery's equation where the ellipsoid shape factor results $\lambda = \kappa$.

If the reaction coefficient $\beta$ in Eq. (35) is large enough $\kappa \approx 0$ implying $\mathbf{p} = \Omega \cdot \mathbf{p}$, that is, tubes rotate according to the flow vorticity, and then the order parameter remains unchanged (the microstructure is rotating but the magnitude of the orientation remains unchanged).

Thus, if we come back to the orientation/aggregation model described in the previous section, more than increasing artificially the diffusion coefficient to resist the orientation induced by the Jeffery's kinematics related to infinite aspect ratio tubes, one could simply consider that in the Jeffery's equation the ellipsoids factor evolves in a way allowing tubes in aggregates to have $\lambda = 0$ allowing the rotation of the microstructure without modifying the order parameter that characterizes the magnitude of the orientation. We propose a model in this direction in the next section.

## A two population model

In this section we introduce a model that consists of two populations of tubes inspired of models proposed for associative polymers [8, 21]. One related to the entangled tubes constituting the aggregates, and the other one that concerns isolated tubes that can freely align according to the Jeffery's equation for particles with infinite aspect ratio. In what follows, we assume a homogeneous system whose microstructure does not depend on the physical space.

The distribution functions related to these two populations are denoted by $\Xi(\mathbf{p})$ and $\Gamma(\mathbf{p})$ for entangled and free tubes respectively. We are assuming that aggregates are being created at rate $v_c$ and destroyed at rate $v_d$, both depending mainly on the flow shear rate. These distribution functions verify the normality condition:

$$ \int_{S(0,1)} (\Xi(\mathbf{p}, t) + \Gamma(\mathbf{p}, t)) d\mathbf{p} = 1, \quad \forall t \quad (41) $$

The Fokker-Planck equations related to both populations writes:

$$ \left\{ \begin{array}{l} \frac{d\Xi(\mathbf{p})}{dt} + \frac{\partial}{\partial \mathbf{p}} \left\{ \Xi(\mathbf{p}) \frac{d\mathbf{p}}{dt} \right\}_A \\ \frac{d\Gamma(\mathbf{p})}{dt} + \frac{\partial}{\partial \mathbf{p}} \left\{ \Gamma(\mathbf{p}) \frac{d\mathbf{p}}{dt} \right\}_F \end{array} \right. \begin{array}{l} \frac{\partial}{\partial \mathbf{p}} \left\{ D_r \frac{\partial \Xi(\mathbf{p})}{\partial \mathbf{p}} \right\} = v_c \Gamma(\mathbf{p}) \quad v_d \Xi(\mathbf{p}) \\ \frac{\partial}{\partial \mathbf{p}} \left\{ D_r \frac{\partial \Gamma(\mathbf{p})}{\partial \mathbf{p}} \right\} = v_c \Gamma(\mathbf{p}) + v_d \Xi(\mathbf{p}) \end{array} \quad (42) $$

fully compatible with Eq. (41) and where for the sake of simplicity the rotary diffusions of both populations are considered the same. In Eq. (42) the rotary velocities read:

$$ \left\{ \begin{array}{l} \frac{d\mathbf{p}}{dt} \Big|_A = \Omega \cdot \mathbf{p} \\ \frac{d\mathbf{p}}{dt} \Big|_F = \nabla \mathbf{v} \cdot \mathbf{p} - (\mathbf{p} \cdot \mathbf{D} \cdot \mathbf{p}) \mathbf{p} \end{array} \right. \quad (43) $$

that agrees with the experimental observations that reveal that CNT aggregates rotate without deforming.

Now, the extra-stress tensor related to the presence of tubes writes:

$$ \tau = 2\eta \left( (N_p)_A ({}^A \mathbf{a}_4 : \mathbf{D}) + (N_p)_F ({}^F \mathbf{a}_4 : \mathbf{D}) \right) \quad (44) $$

where $(N_p)_A$ and $(N_p)_F$ are two rheological parameters related to the contribution of aggregates and free tubes to the stress respectively. The two fourth order orientation tensors in Eq. (45) are given by:

$$ \left\{ \begin{array}{l} ^A \mathbf{a}_4 = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \Xi(\mathbf{p}) d\mathbf{p} \\ ^F \mathbf{a}_4 = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \Gamma(\mathbf{p}) d\mathbf{p} \end{array} \right. \quad (45) $$

If we define the fraction of both populations of tubes:

$$ \left\{ \begin{array}{l} C^A = \int_{S(0,1)} \Xi(\mathbf{p}) d\mathbf{p} \\ C^F = \int_{S(0,1)} \Gamma(\mathbf{p}) d\mathbf{p} \end{array} \right. \quad (46) $$

the integration on the conformation space of Eq. (42) results:

$$ \left\{ \begin{array}{l} \frac{dC^A}{dt} = v_c C^F - v_d C^A \\ \frac{dC^F}{dt} = -v_c C^F + v_d C^A \end{array} \right. \quad (47) $$

Adding both expressions in Eq. (47) we obtain the conservation balance:

$$ \frac{dC^A}{dt} + \frac{dC^F}{dt} = 0 \quad (48) $$

On the other hand because the normality condition:

$$ C^A + C^F = 1 \quad (49) $$
Now, in the steady state, Eqs. (47) and (49) lead to the equilibrium fractions:

$$ \left\{ \begin{aligned} C^A &= \frac{v_c}{v_c+v_d} \\ C^F &= \frac{v_d}{v_c+v_d} \end{aligned} \right. \qquad (50) $$

that reduces as expected to $C^A \approx 1$ ($C^F \approx 0$) if $v_c \gg v_d$; and $C^F \approx 1$ ($C^A \approx 0$) if $v_d \gg v_c$.

One could expect that $(N_p)_A \gg (N_p)_F$ (the contribution of aggregates to the stress is much higher than the one due to the free tubes, and the aggregates are always present in processes) from which Eq. (44) reduces to:

$$ \tau \approx 2\eta(N_p)_A(^A\mathbf{a}_4 : \mathbf{D}) \qquad (51) $$

If a simple shear flow with a very low shear rate is applied to the suspension ($v^T = (\Gamma_y, 0, 0), \Gamma \gg 1$) we can

assume that the resulting steady state will be characterized by: $C^A \approx 1$ and $\Xi \approx \frac{1}{4\pi}$. Thus, if we define $a = (^A\mathbf{a}_{1212}) = \int_{S(0,1)} \mathbf{p}_1^2 \cdot \mathbf{p}_2^2 \frac{1}{4\pi} d\mathbf{p}$, the shear stress writes:

$$ \tau_{12} = 2\eta\left(\frac{1}{2} + (N_p)_A a\right)\Gamma \qquad (52) $$

which allows computing the rheological parameter $(N_p)_A$ from the plateau of the flow curves.

The remaining rheological parameters, the aggregation and disaggregation velocities, can be identified by fitting the flow curves at moderate shear rates. For this purpose the global model should be solved. Two modeling possibilities exist: (i) a fully kinetic theory modeling, accurate but computationally expensive, or (ii) a macroscopic model making use of the orientation tensors.

<figure><img src="image_6.png" /><figcaption>Fig. 5 Orientation field in the $r \theta$ plane. Time steps 1, 10, 20 and 30. $\Delta t = 0.01$s. $\lambda = 1$</figcaption></figure>
Kinetic theory model

The resulting kinetic theory model is defined by:

$$ \left\{ \begin{array}{l} \frac{d\Xi(\mathbf{p})}{dt} + \frac{\partial}{\partial \mathbf{p}} \left\{ \Xi(\mathbf{p}) \frac{d\mathbf{p}}{dt} \right\}_A = v_c \Gamma(\mathbf{p}) \quad v_d \Xi(\mathbf{p}) \\ \frac{d\Gamma(\mathbf{p})}{dt} + \frac{\partial}{\partial \mathbf{p}} \left\{ \Gamma(\mathbf{p}) \frac{d\mathbf{p}}{dt} \right\}_F = v_c \Gamma(\mathbf{p}) + v_d \Xi(\mathbf{p}) \\ {}^A\mathbf{a}_4 = \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \Xi(\mathbf{p}) d\mathbf{p} \\ \tau \approx 2\eta(N_p) {}^A\mathbf{a}_4 : \mathbf{D} \end{array} \right. \quad (53) $$

Macroscopic modelling

To avoid the solution of the Fokker-Planck equations we could derive a macroscopic model making use of different

orientation tensors: (i) the fourth order ones defined by Eq. (45) and the second order ones given by

$$ \left\{ \begin{aligned} {}^A\mathbf{a}_2 &= \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \Xi(\mathbf{p}) \, d\mathbf{p} \\ {}^F\mathbf{a}_2 &= \int_{S(0,1)} \mathbf{p} \otimes \mathbf{p} \Gamma(\mathbf{p}) \, d\mathbf{p} \end{aligned} \right. \quad (54) $$

We could derive the evolution equation of both second order orientation tensors by taking the time derivative of both expressions in Eq. (54) and then introducing Eqs. (43) and (42).

First, we consider:

$$ \frac{d^A \mathbf{a}_2}{dt} = \int_{S(0,1)} \left( \dot{\mathbf{p}} \otimes \mathbf{p} \Xi(\mathbf{p}) + \mathbf{p} \otimes \dot{\mathbf{p}} \Xi(\mathbf{p}) + \mathbf{p} \otimes \mathbf{p} \frac{d\Xi(\mathbf{p})}{dt} \right) d\mathbf{p} \quad (55) $$

<figure><img src="image_7.png" /><figcaption>Fig. 6 Orientation field in the r θ plane. Time steps 1, 10, 20 and 30, with Δt = 0.01s. λ = 0. Note that in this case the orientation field remains more isotropic than in the preceding case</figcaption></figure>
that leads to:

$$ \frac{d^A \mathbf{a}_2}{dt} = \mathbf{\Omega} \cdot \mathbf{a}_2 - \mathbf{a}_2 \cdot \mathbf{\Omega} - 6D_r \left( \mathbf{a}_2 - \frac{\mathbf{I}}{3} \right) + v_c^F \mathbf{a}_2 - v_d^A \mathbf{a}_2 \quad (56) $$

that does not involve any closure relation. The evolution of the second order orientation tensor related to the free tubes population reads:

$$ \frac{d^F \mathbf{a}_2}{dt} = \nabla \mathbf{v} \cdot \mathbf{F} \mathbf{a}_2 - \mathbf{F} \mathbf{a}_2 \cdot (\nabla \mathbf{v})^T - 2 \mathbf{F} \mathbf{a}_4 : \mathbf{D} - 6D_r \left( \mathbf{F} \mathbf{a}_2 - \frac{\mathbf{I}}{3} \right) - v_c^F \mathbf{a}_2 + v_d^A \mathbf{a}_2 \quad (57) $$

whose solution needs the introduction of a closure relation, being the simplest ones the quadratic closure $\mathbf{F} \mathbf{a}_4 = \mathbf{F} \mathbf{a}_2 \otimes \mathbf{F} \mathbf{a}_2$, the linear or the hybrid ones. In the numerical examples presented in the next section the quadratic closure relation was considered.

The resulting model consists of Eqs. (56) and (57) combined with an appropriate closure relation and the constitutive Eq. (44).

## Numerical experiments

In this section we are simulating spin coating processes. For this purpose we are coupling the flow kinematics with the macroscopic modelling just described involving two populations and the evolution of the associated orientation tensors. We are simply illustrating the model solution, the model calibration using rheometric flows is a work in progress.

The numerical technique employed for solving the momentum balance equations is the Natural Element Method (NEM), deeply described in Gonzalez et al. [7]. In Cueto et al. [4], we described the application of such a technique for simulating the free surface problem related to spin coating processes.

In what follows we are considering the two limit cases: (i) the one in which the aggregation velocity is much higher than the disaggregation one and (ii) the opposite case in which the disaggregation velocity is much higher than the aggregation one. For all the cases, the before presented model was applied to the simulation of spinning a drop of CNT suspension, with the assumption that the initial shape is Gaussian. The model was composed of 977 nodes under axisymmetric assumptions for the flow (but not for the orientation field, which is three-dimensional). Non-slip boundary conditions were assumed at the horizontal plane, whereas symmetry boundary conditions are assumed at the axis of symmetry. The initial orientation distribution was in all cases isotropic.

Disaggregation velocity bigger than aggregation one

If we consider that the disaggregation velocity is always dominant over the aggregation one, the resulting behaviour of the fluid will be essentially the same described in Cueto et al. [4]. Here we consider the orientation as governed by a factor $\lambda = 1$. The results concerning the CNTs orientation are shown in Fig. 5a–d). It can be noticed how the CNTs tend to rapidly align with the flow. To understand the apparent isotropic final distribution of the CNTs it must be noticed that the orientation tensor is fully three-dimensional, and that CNTs tend to orient with the divergent flow field, and thus in the r-θ plane. That is why the ellipses are of considerable lower radius than those at initial time steps.

Aggregation velocity bigger than the disaggregation one

In this case, if we consider that the aggregation velocity is dominant, aggregates dominate the resulting behaviour of the suspension, and therefore tend to avoid the orientation dictated by the flow kinematics. We have simulated the same spin coating process with a parameter $\lambda = 0$. Results are depicted in Fig. 6a–d. The computed orientation state is in agreement with our expectations.

## Conclusions

In this paper a new model for CNT suspensions is developed. It deals with non-treated CNTs, and therefore with the possibility of aggregation of the nanotubes. If aggregates are formed the proposed model predicts a lower tendency to alignment within the imposed flow, if compared to the one dictated by Jeffery's equation. This has been observed experimentally also.

From the numerical point of view, the model has been solved by employing Natural Element methods, a technique that is able to easily track free surface flows. An application to spin coating processes has been developed. It has been noticed how the bigger the presence of aggregates is, the lower tendency to align to the flow is observed in the simulations.

**Acknowledgements** This work has been partially supported by the Spanish Ministry of Science and Innovation, through grant number CICYT DPI2008 00918. This support is gratefully acknowledged.

## References

1. Advani SG, Tucker CL III (1990) Closure approximations for three dimensional structure tensors. J Rheol 34:367 386

2. Ammar A, Mokdad B, Chinesta F, Keunings R (2006) A new family of solvers for some classes of multidimensional partial
differential equations encountered in kinetic theory modelling of
complex fluids. J Non Newton Fluid Mech 139:153 176

3. Chinesta F, Chaidron G, Poitou A (2003) On the solution of the Fokker Planck equation in steady recirculating flows involving short fibre suspensions. J Non Newton Fluid Mech 113:97 125

4. Cueto E, Ma A, Chinesta F, Mackley M (2008) Numerical simulation of spin coating processes with functionalized carbon nanotubes suspensions. Int J Mater Form 2:89 99

5. Dupret F, Verleye V (1999) Modelling the flow of fibre suspensions in narrows gaps. In: Siginer DA, De Kee D, Chabra RP (eds) Advances in the flow and rheology of non Newtonian fluids. Rheology series. Elsevier, Amsterdam, pp 1347 1398

6. Folgar F, Tucker CL III (1984) Orientation behaviour of fibers in concentrated suspensions. J Reinf Plast Compos 3:98 119

7. Gonzalez D, Cueto E, Chinesta F, Doblare M (2007) An updated Lagrangian strategy for free surface fluid dynamics. J Comput Phys 223:127 150

8. Hernandez Cifre JG, ThMAOM Barenbrug, Schieber JD, van den Brule BHAA (2003) Brownian dynamics simulation of reversible polymer networks under shear using a non interacting dumbbell model. J Non Newton Fluid Mech 113:73 96

9. Hobbie EK, Fry DJ (2007) Rheology of concentrated carbon nanotube suspensions. J Chem Phys 126:124907 1 124907 7

10. Jeffery GB (1922) The motion of ellipsoidal particles immersed in a viscous fluid. Proc R Soc London A102:161 179

11. Kinloch IA, Roberts SA, Windle AH (2002) A rheological study of concentrated aqueous nanotube dispersions. Polymer 43:7483 7491

12. Lin Gibson S, Pathak JA, Grulke EA, Wang H, Hobbie EK (2004) Elastic flow instability in nanotube suspensions. Phys Rev Lett 92:048302 1 048302 4

13. Ma AWK, Mackley MR, Rahatekar SS (2007) Experimental observation on the flow induced assembly of carbon nanotube suspensions to form helical bands. Rheol Acta 46:979 987

14. Ma AWK, Chinesta F, Mackley MR (2008a) Rheological modelling of carbon nanotube aggregate suspensions. J Rheol 52(6):1311 1330

15. Ma AWK, Chinesta F, Tuladhar T, Mackley MR (2008b) Filament stretching of carbon nanotube suspensions. Rheol Acta 47:447 457

16. Ma AWK, Mackley MR, Chinesta F (2008c) The microstructure and rheology of carbon nanotube suspensions. Int J Mater Form, published online

17. Ma AWK, Chinesta F, Mackley MR (2009) The rheology and modelling of chemically treated carbon nanotube suspensions. J Rheol 53(3):547 573

18. McNally T, Pöschke P, Halley P, Murphy M, Martin D, Bell SEJ, Brennan GP, Bein D, Lemoine P, Quinn JP (2005) Polyethylene multiwalled carbon nanotube composites. Polymer 46:8222 8232

19. Petrie CJS (1999) The rheology of fibre suspensions. J Non Newton Fluid Mech 87:369 402

20. Rahatekar SS, Koziol KKK, Butler SA, Elliott JA, Shaffer MSP, Mackley MR, Windle AH (2006) Viscosity enhancement of an epoxy resin matrix by addition of multi wall carbon nanotubes. J Rheol 50:599 610

21. Vaccaro A, Marrucci G (2000) A model for the nonlinear rheology of associating polymers. J Non Newton Fluid Mech 92:261 273

22. Xu J, Chatterjee S, Koelling KW, Wang Y, Bechtel SE (2005) Shear and extensional rheology of carbon nanofiber suspensions. Rheol Acta 44:537 562