The rheology and modeling of chemically treated carbon nanotubes suspensions

his paper reports recent experimental findings and rheological modeling on chemically treated ingle-walled carbon nanotubes CNTs suspended within an epoxy resin. When a CNT suspension as subject to a steady shear flow, it exhibited a shear-thinning characteristic, which was ubsequently modeled by a Fokker–Planck FP based orientation model. The model assumes that he shear flow aligns CNT in the flow direction, but there are events such as Brownian motion and ube–tube interaction trying to randomize the orientation. In the FP orientation model, randomizing vents were modeled with an appropriate rotary diffusion coefficient Dr and the shear-thinning ehavior was explained in terms of progressive alignment of CNTs toward the shear direction. In erms of linear viscoelasticity LVE , small-amplitude oscillatory measurements revealed mild lasticity for semidilute treated CNT suspensions. The exact origin for this elasticity is not clear nd both tube–tube interaction and bending/stretching of CNTs have been proposed by other uthors as possible origins. It is, however, clear from the current modeling that the experimental volution of storage modulus G cannot be described using a single-mode Maxwell model or imple Brownian rod modeling. In this paper, experimental LVE data of the treated CNT uspensions were fitted using the FP orientation model with an “effective diffusion coefficient” erm and an empirical relation was subsequently identified for the effective diffusion term. ntuitively, chemical treatment has created a weakly interconnected network of CNT and it is elieved that the mild elasticity originated from this weak network as well as other randomizing vents Brownian motion and tube–tube hydrodynamic interaction . Finally, step strain xperiments confirmed the presence of a weak network at small strains, which at large strains was ound to be destroyed. Incorporation of a strain softening factor allowed for the formulation of a


I. INTRODUCTION
Carbon nanotubes ͑CNTs͒ are cylinders of rolled graphite sheets possessing high mechanical strength, low density, and special electronic properties ͓Iijima ͑1991͒; Dresselhaus and Dai ͑2004͔͒.They belong to a relatively new class of fibrous material that can potentially be used for high-performance nanocomposites Huang et al. ͑2006͒ reported that a mixing time of several hours was required to establish a stable rheology for CNT suspensions and even then uniformity at nanoscale is not guaranteed.In this regard, introducing functional groups onto the wall/cap of CNTs also known as functionalization has proven to be one of the most effective ways to enhance dispersion of CNT in a processing matrix ͓Kinloch et al. ͑2002͒; Banerjee et al. ͑2003͒; Dyke and Tour ͑2003͒; Dyke and Tour ͑2004͔͒.Despite the potential of functionalization, a limited amount of work has been carried out on the rheological modeling of chemically treated CNT suspensions.This paper presents experimental results and modeling of chemically treated CNTs suspended in an epoxy resin, considering both steady shear and linear viscoelastic ͑LVE͒ responses for this type of suspension.

II. EXPERIMENTAL DETAILS
CNTs used were single-walled CNTs produced by high pressure carbon monoxide disproportionation ͑HiPco͒ method and they were supplied by Nanocomposites Inc., USA.The exact length and diameter distribution for the treated CNTs have not been fully characterized; however, according to the supplier, each individual CNT was 0.5-1 m long and with a typical diameter of 1.2 nm.In the case of treated CNTs, effectivestacking between CNTs and subsequent aggregation was prevented by covalently attaching arene diazonium salts onto the sidewall of CNTs ͓Dyke and Tour ͑2003͒; Dyke and Tour ͑2004͔͒.For comparison, some examples of basic rheology and microstructure of untreated CNTs are also considered in this paper.0.6 g of CNTs was dispersed in 60 g of epoxy resin ͑Araldite LY556, Huntsman Inc.͒ to give a 1% masterbatch suspension.Suspensions with lower CNT concentrations were obtained by a dilution of the masterbatch.Mixing was carried out for 5 h using a Silverson-L4R homogenizer and the microstructure of resulting mixtures was optically characterized using the Cambridge shear system ͑Linkam Scientific Instruments Ltd.͒ ͓Bower et al. ͑1998͒; Mackley et al. ͑1999͔͒.Figure 1 compares the typical optical texture of a suspension with chemically treated CNTs and that with untreated CNTs.Some impurities of amorphous carbon and metal catalyst residue were present in the chemically treated CNT suspension ͓Fig.1͑a͔͒; however it was clear from the micrograph that the suspension showed no optically resolvable aggregates of CNTs, and the mixture was well dispersed at the micronlevel.In addition, no formation of optically resolvable aggregates was observed in the presence of shear flow, as shown in Fig. 2. By contrast, the untreated CNT suspension consisted of optically resolvable CNT aggregates ͓Fig.1͑b͔͒ as discussed by Rahatekar et al. ͑2006͒ andMa et al. ͑2008a͒.
The extensional rheology of the same treated CNT sample was characterized and reported by Ma et al. ͑2008c͒.By fitting the experimental extensional viscosity data to the models of Batchelor ͑1971͒ and Shaqfeh and Fredricksen ͑1990͒, the average aspect ratio was estimated to be 180, which is in reasonable agreement with the atomic force microscopy ͑AFM͒ results where small bundles with a typical diameter of 5 nm were observed after high shear mixing.It should, however, be noted that the CNTs were not "monodispersed" and there was a certain length ͑and diameter͒ distribution associated with the CNTs used in the actual experiments-originated from the synthesis, mixing, and subsequent chemical treatment processes ͓Ma et al. ͑2008b͔͒.In this study, the aspect ratio distribution of CNT has not been studied in detail.In terms of modeling, given an average aspect ratio in the order of 180, the shape factor k in the Jeffery equation ͓Eq.͑4͔͒ was assumed to be 1.The aspect ratio distribution could have an influence on the overall rheological behavior of the system ͑as discussed, for example, in the paper of Pittman and Bayram ͑1990͔͒, but the effect of length distribution has not been considered explicitly in the current modeling.If there is any effect due to the "polydispersity" of the CNTs, this would have been accounted for in the subsequent model fitting, as represented by the parameter N p in Eq. ͑12͒, which combines the effect of concentration and aspect ratio.Rheological measurements were made using an Advanced Rheometric Expansion System ͑ARES͒ strain-controlled rheometer with 50 mm parallel plates and a gap size of 0.3 mm.The temperature was maintained at 25 °C by an air oven.In steady shear experiments, the apparent viscosity of the sample was measured after a particular shear rate was applied for 100 s.In the small-amplitude oscillatory shear experiment, a strain amplitude of 1% was used.In order to minimize possible complication from sample loading, samples were squeezed between the parallel plates slowly and were rested for at least 2 h before any measurements were carried out.
Step strain experiments were carried out in order to explore the transition from small to large strain deformations.

A. Experimental results
Steady shear experimental plots of the apparent viscosity ͑ a ͒ as a function of shear rate ͑␥ ˙͒ for chemically treated and untreated CNT suspensions are shown in Fig. 3.Both suspensions were shear thinned asymptotically to the matrix viscosity ͑also shown͒, but given the same weight concentration of 0.33%, the suspension with untreated CNTs showed a significantly larger viscosity enhancement effect than the suspension with chemically treated CNTs.Previous work ͓Rahatekar et al. ͑2006͔͒ identified the viscosity enhancement for untreated CNT suspensions with the optically observed aggregate microstructure and similar observations were also reported in a comparative study of chemically treated and untreated carbon nanofiber ͑CNF͒ suspensions ͓Xu et al. ͑2005͔͒.The data obtained in Fig. 3 were found to be consistent irrespective of whether a particular shear rate was approached from either a higher or lower shear rate provided an adequate equilibration time was allowed at each respective shear rate.

B. Fokker-Planck based simple-orientation model
In order to model the steady shear-thinning response of treated CNTs, an orientation model was developed.In this model, the CNTs were assumed to be an assembly of high aspect ratio rigid fibers subjected to shear and rotary diffusion.Similar models have been used to model short fiber suspensions; however, the numerical values of key parameters in the current modeling were identified by fitting to experimental data and no theoretical relations as described, for example, in Petrie ͑1999͒, were used.The model proposed is consistent with the shear-induced alignment of CNTs as observed using optical microscopy ͓Rahatekar et al. 2006͔͒ and transmission electron microscopy ͓Fan and Advani ͑2005, 2007͔͒.In terms of the mathematical treatment of CNT orientation distribution, a Fokker-Planck ͑FP͒ description, instead of the readily available closure approximations for the orientation tensors, was used in the steady shear modeling.Nonlocal effects such as tube-tube interaction and CNT-polymer interaction can be introduced as proposed in the work of Forest and Wang ͑2005͒, where they considered blends of rodlike liquid crystal polymer and flexible polymers.The latter approach, however, requires some additional information that is not readily available.This includes detailed descriptions for ͑i͒ tube-tube interaction, ͑ii͒ CNT conformation, and ͑iii͒ CNT-polymer interaction.Besides that, taking into consideration nonlocal interactions would lead to extraelastic body forces that cannot be written as the divergence of a second order tensor and therefore cannot be directly introduced in the general expression for the stress tensor.As pointed out by Forest and Wang ͑2005͒, it is necessary to introduce the extraelastic body forces in the momentum balance equation, making this type of analysis more complicated to perform.This study aims to derive a simple constitutive model for chemically treated CNT suspensions and for this reason, tube-tube interaction is represented using an effective diffusion coefficient and CNT-polymer interaction is neglected in our current system where the suspending matrix is made up of small prepolymer molecules.Experimentally, similar steady shear and LVE rheological responses were obtained when the epoxy matrix was replaced by an acrylic resin and this further confirms the belief that CNT-polymer interaction plays a relatively weak role in terms of the overall rheology of the chemically treated CNT suspensions.

Model formulation
The momentum balance equation neglecting the mass and inertia terms reads as div͑ =͒ = 0, ͑1͒ where = is the total stress tensor.The mass balance equation for an incompressible fluid gives div͑v គ ͑x គ͒͒ = 0, ͑2͒ where v គ ͑x គ͒ is the velocity field.
The constitutive equation for a dilute suspension with high aspect ratio rigid fibers can be written as ͓Batchelor ͑1970͒; Hinch and Leal ͑1975, 1976͒; Petrie ͑1999͔͒ where P denotes the pressure, I = is the unit tensor, is the suspending medium viscosity, D = is the strain rate tensor ͑symmetric component of the velocity gradient tensor͒, and = f is the stress contribution due to the presence of fibers.
The suspension microstructure can be defined using the fiber orientation.It was assumed that the fibers are ellipsoidal and rigid and that they are immersed in a Newtonian matrix whose kinematics are defined by grad͑v គ ͑x គ͒͒.If p denotes the unit vector aligned in the fiber axis direction, then, its evolution will be given by the Jeffery equation ͓Jeffery ͑1922͔͒, where k = ͑ 2 −1͒ / ͑ 2 +1͒, is the ͑effective͒ fiber aspect ratio, and ⍀ = denotes the vorticity tensor.The tensor product of vectors a and b is defined as ͑a គ b គ ͒ ij = a i b j , and ":" denotes the tensor product contracted twice ͓i.e., ͑a = : b =͒ = ͚͚a ij b ji ͔.
Although the Jeffery equation was first proposed for an isolated fiber immersed in a Newtonian fluid, it can also be applied to a population of fibers if the suspension is dilute and the interaction between the fibers is negligible.However, defining the fluid microstructure using the orientation of each fiber in the suspension is not helpful for the formulation of a mesoscale model.Instead, a more useful way of describing the microstructure is to use a kinetic theory approach which introduces an orientation distribution function ͑x គ , p គ , t͒ such that ͑x គ , p គ , t͒dp គ represents the probability of finding at point x គ and time t: a fiber whose orientation is within the interval defined by p គ and p គ + dp គ .In this expression, the physical coordinates ͑x គ , t͒ ͑space and time͒ can be distinguished from the conformation coordinate p គ , which describes the orientation defined on the surface of the unit sphere S͑0,1͒.This distribution function satisfies the normality condition, The consideration of an orientation distribution function allows certain moments to be defined: namely, the second order moment, also known as the second order orientation tensor, and the fourth order moment also known as the fourth order orientation tensor, The evolution of the orientation distribution function is governed by the FP equation ͓see, for example, Hinch and Leal ͑1972͔͒,

͑8͒
where the advection field dp គ / dt is given by the Jeffery's equation ͓Eq.͑4͔͒, and d / dt is the material derivative, i.e., D r given in Eq. ͑8͒ is the rotary diffusion coefficient of the fibers and is a key adjustable parameter within the model to describe both orientation and randomizing events in semidilute suspensions.Folgar and Tucker ͑1984͒ proposed the use of an empirical effective rotary diffusion ͑D ˆr͒ and in their treatment, D ˆr is proportional to the shear rate.
As the flow kinematics induce a homogeneous distribution state, the material derivative reduces to the partial derivative, which vanishes when looking for the steady solution.In the present work, it is assumed that the orientation does not depend on the physical coordinates, and therefore = ͑p គ , t͒.By integrating Eq. ͑8͒ in the conformation space ͓the unit sphere surface S͑0,1͔͒, it is easy to verify that

͑11͒
which implies that if Eq. ͑8͒ is solved using an initial distribution that satisfies the normality condition, its time evolution will also verify the normality condition.On the other hand, if one is looking directly for the steady orientation distribution, the normality condition related to the searched distribution should be imposed using an appropriate technique ͑e.g., Lagrange multipliers͒.
An appropriate way to derive a constitutive equation for the suspension involves the use of spatial homogenization and statistical averaging.Together with other simplifying hypotheses for high aspect ratio fibers, the stress tensor = f due to the presence of fibers can be defined ͓Batchelor ͑1971͒; Hinch and Leal ͑1975͒; Hinch and Leal ͑1976͔͒, where N p is a scalar parameter that depends on the fiber concentration as well as the aspect ratio of fibers in dilute suspensions whereas for dilute suspensions ␤ =6 str ͑ is the number concentration and str is the viscous drag coefficient͒ ͓Larson ͑1999͔͒.The second term on the right-hand side of Eq. ͑12͒ essentially accounts for the effect of Brownian motion; but in general, other randomizing events, such as fiber-fiber hydrodynamic interaction, can be accounted for by replacing D r with an "effective diffusion term."In terms of steady shear or large strain deformations, the second term becomes negligible compared with the first term involving the fourth order orientation tensor.The second term in Eq. ͑12͒ was, however, found to be important in explaining the mild elasticity observed in small-strain LVE measurements.The second term is included in the LVE modeling, but it can be neglected in steady shear data fitting.It should be noted that although the expression proposed by Hinch and Leal ͑1975, 1976͒ is used in the current study, there exist alternative formulae as reported by Chan and Terentjev ͑2007͒ and Grmela ͑2008͒.
Instead of using a kinetic theory model ͓Eqs.͑3͒, ͑7͒, ͑8͒, and ͑12͔͒, some authors prefer to use purely macroscopic models because of lower computation requirements.The simplest macroscopic model can be derived by taking the time derivative of Eq. ͑6͒ and then introducing both the Jeffery equation ͓Eq.͑4͔͒ and the FP equation ͓Eq.͑8͔͒.Prager ͑1957͒ showed that As shown in the equation above, the evolution of the second order orientation tensor ͑a =͒ depends on the fourth order orientation tensor ͑a = = ͒ and closure relations that express the fourth order orientation tensor as a function of the lower order orientation tensors is commonly used to solve the mathematical problem in a closed form.Different forms of closure relationship, including quadratic, linear, hybrid and natural relations, can be found in the literature ͓Advani and Tucker ͑1990͒; Dupret and Verleye ͑1999͔͒.However, it is worth noting that these closure relations have their own constraints and the use of these relations could sometimes adversely influence the computed numerical solutions.
For instance, the quadratic closure relation becomes exact only for perfectly aligned fibers whereas the linear closure relation becomes exact only for randomly orientated fibers ͓Petrie ͑1999͔͒.In steady shear modeling, the orientation distribution of CNT was expected to be neither completely isotropic nor randomly orientated ͑for most shear conditions͒; closure approximations were therefore avoided and the FP equation was solved directly instead ͓Öttinger and Laso ͑1992͒; Once the distribution function is known, the fourth order orientation tensor and the stress tensor can then be computed.For simple shear flow and with an assumption that the flow is not perturbed by the presence and orientation of fibers,

͑14͒
In the case of steady shear or large strain deformation, the stress contribution from Brownian motion and other randomizing effects is negligible compared with the anisotropic viscous term involving the fourth order orientation tensor; the shear stress according to Eq. ͑12͒ is therefore given by 12 = ␥ ˙+ 2N p a 1212 ␥ ˙͑15͒ and the apparent viscosity is therefore

͑16͒
The described model was originally developed for dilute suspensions, it has however been used in several commercial simulation codes and has been successfully applied to model semi-dilute and concentrated suspensions generally encountered in industrial applications ͓Advani ͑1994͔͒.In such cases, the model previously described could still be used to compute the flow kinematics and the microstructure evolution, but rheological parameters in the model, such as N p and D r , should be determined from experimental data.As an alternative, there also exists models that allows for the estimation of D r and N p ͓see, for example, Larson ͑1999͔͒, but the hypothesis of which they are based on is only valid under certain circumstances.

C. Model fitting
The FP based simple-orientation model has been applied to the treated CNT suspensions and Fig. 4͑a͒ shows the model fitting of chemically treated CNT suspensions with four different weight concentrations ͑0.05%, 0.2%, 0.33%, and 0.5%͒.Taking the 0.33% CNT suspension as an example, the evolution of a in the orientation model is controlled by the values of D r and N p and the best constant-D r fit was obtained for D r = 0.005 s −1 and N p = 7.In Fig. 4͑b͒, different values of D r were used to illustrate the sensitivity of the model on the fitting parameter D r .In general, there is a good agreement between a predicted by the model and the experimental data; but for high shear rates ͑␥ ˙͒, the ͓Note: by setting N p = 0, the fiber contribution term in Eq. ͑16͒ becomes zero and therefore, a = where is the suspending medium viscosity.͔predicted a is slightly higher than the experimental value.For instance, at ␥ ˙=60 s −1 , a of the 0.33% suspension was found experimentally to be 10.5 Pa s, but the orientation model predicted a viscosity of 11.3 Pa s ͑with a constant D r = 0.005 s −1 ͒ and this gives an error of about 8%.Although there is an error of a few percents in predicting high shear-rate data, the use of a constant D r should be sufficient to provide a reasonably good estimation of a for general engineering problems, and this also saves the effort in determining the exact relationship between D r and ␥ ˙.However, in cases where a more accurate description of high shear viscosity is needed, the dependence of D r on ␥ ˙should be carefully evaluated ͓see, for example, Folgar and Tucker ͑1984͒ and Larson ͑1999͔͒.
In terms of low shear-rate data, although the torque generated by ␥ ˙Ͻ 0.1 s −1 was too small to be measured accurately and corresponding a data was not available experimentally, the model predicted a low shear-rate Newtonian plateau and the level of the plateau depends on the parameter N p ͓Fig.4͑a͔͒.As shown in the inset figure, N p scaled linearly with the CNT concentration.
Given the best-fit value for D r to be 0.005 s −1 , one can compare this value to the theoretical prediction by Doi and Edwards ͑1986͒ for dilute suspensions where interparticle interactions are absent and rotary diffusion is due to Brownian motion only, where k B is the Boltzmann constant ͑1.38ϫ 10 −23 J / K͒, T is the temperature ͑in Kelvin͒, L is the CNT length, d is the CNT diameter, and is the viscosity of the suspending medium.
For the treated CNTs having an aspect ratio ͑L / d͒ of about 180 and epoxy resin with a base viscosity of 10 Pa s, Eq. ͑17͒ predicted a Brownian rotary diffusion ͑D ro ͒ of 0.002 s −1 , which is in the same order of magnitude as the D r determined from the experimental data fitting.The higher value of D r could probably be explained by extrainterparticle interactions in the case of semidilute suspensions.
To quantify the degree of CNT alignment, a scalar order parameter ͑S͒ can be used ͓see, for example, Fry et al. ͑2005͔͒ and the definition can be found in many liquid crystalline polymer textbooks such as Donald et al. ͑2006͒ and Laso et al. ͑2006͒, where Q = = ͑a = − I = / 3͒, a = is the second order orientation tensor as defined in Eq. ͑6͒, and I = is the identity matrix.S is a scalar parameter varies between 0 and 1 and the larger the magnitude of S, the higher the degree of alignment.Figure 4͑c͒ presents the scalar order parameter ͑S͒ as a function of shear rate as predicted by the orientation modeling.At a shear rate of 100 s −1 , the model predicted an order parameter of 0.92.
The orientation model, in principle, could be used to model CNT suspensions where the CNTs do not aggregate in the presence of shear flow.However, recent modeling work carried out by us ͓Ma et al. ͑2008a͔͒ suggests that such a model cannot satisfactorily describe the shear-thinning characteristics for ͑untreated͒ CNT suspensions where optically resolvable CNT aggregates are present ͓see, for example, Rahatekar et al. ͑2006͔͒.For CNT suspensions where the CNTs aggregated in shear flow, a more advanced model named the "aggregation/orientation ͑AO͒ model" has been developed to explain the shear-thinning behavior observed experimentally.The AO model considered a hierarchy of aggregate structures within the CNT suspensions and the evolution of CNT orientation depended on the flow conditions as well as the entanglement state of CNTs ͓Ma et al. ͑2008a͔͒.

A. Experimental results
LVE of CNT suspensions was studied using small-amplitude oscillatory measurements.Epoxy resin showed scattered GЈ data with torque values very close to the detection limit of the transducer, implying that the elasticity of the matrix ͑G epoxy Ј ϳ 0͒ is negligible.Epoxy behaved essentially as a Newtonian fluid with viscous dissipation that is consistent with steady shear measurements.Addition of CNTs increased the values of both GЈ and GЉ as shown in Figs.5͑a͒ and 5͑b͒.Measurements were made at a strain of 1%, which was well within the linear strain response of the suspensions.The enhancement of GЈ was concentration dependent and was more pronounced at high concentration levels ͑0.2% and 0.5%͒.The evolution of GЈ as a function of frequency is consistent with experimental results reported by others ͓Song and Youn ͑2005͒; Xu et al. ͑2005͔͒ and similar evolution of GЈ was also observed experimentally when the treated CNTs were suspended in a Newtonian acrylic resin having a different base viscosity ͓dipropylene glycol diacrylate ͑DPGDA͒; Cytec Industries Inc.͔. Figure 5͑a͒ clearly illustrates that the addition of CNTs increased the elasticity of the system as a whole.This response is very different from that of a typical short fiber suspension, where in the latter case, the addition of non-Brownian fibers was reported to have no extra contribution to the storage modulus ͑GЈ͒ value of the suspending medium ͓Carter ͑1967͒; Ganani and Powell ͑1986͔͒.
To assess the relative importance of viscous and elastic contributions at a given concentration, full GЈ, GЉ, and ‫ء‬ data of the 0.5% CNT suspension are presented in Fig. 5͑c͒ . The value of GЉ was observed to be higher than GЈ within the full range of frequency studied, which implies that the elasticity associated with the addition of CNTs is relatively mild.Although the elastic response is relatively weak, it is interesting to note that the experimental evolution of GЈ and GЉ does not follow the prediction of a singlemode Maxwell model as predicted for Brownian rod systems ͓Larson ͑1999͔͒.

B. Model fitting
In the case of CNFs which are analogous to CNTs, Xu et al. ͑2005͒ carried out experiments and applied elastic and rigid dumbbell models to describe both shear and viscoelastic behavior for treated and untreated CNF suspensions.For treated CNF suspensions, they claimed that a Hookean elastic dumbbell model can successfully capture the LVE response.For the system studied in this paper, the experimental shear thinning was found and has successfully been modeled using the FP orientation diffusion model using a constant diffusion coefficient.To model the LVE of treated CNT suspensions, Eqs.͑12͒ and ͑13͒ are considered in the context of small-amplitude oscillatory measurements.Prior to any oscillatory measurements, CNTs are assumed to be isotropically orientated ͑given that no preshear is applied to the sample͒.In the presence of weak small-amplitude oscillatory flows, the orientation of CNTs would remain very close to an isotropic distribution.This allows for the application of a linear closure relation for the fourth order orientation tensor, which takes the following form ͓see, for example, Advani and Tucker ͑1990͔͒: Given that the small-amplitude oscillatory flow does not perturb significantly the isotropic orientation distribution state, the following approximation holds: Combining Eqs.͑13͒ and ͑20͒ gives which can be rearranged as da dt where a ϵ a = 12 , =1/ 6D r , and =1/ 30D r ϵ / 5.
Equation ͑22͒ takes the form of a Maxwell viscoelastic model and for a smallamplitude oscillatory flow where ␥ ˙͑t͒ = i␥ 0 e it , the real and imaginary components of a ͑i.e., a The shear stress according to Eq. ͑12͒ is and the storage modulus ͑GЈ͒ and the loss modulus ͑GЉ͒ are therefore This analysis is in perfect agreement with the expressions given in Larson ͑1999͒ for high aspect ratio Brownian rods.It is clear from this theoretical analysis that the Brownian motion of particles would result in Maxwell type elasticity, where GЈ ϰ 2 at low frequencies and there is a GЈ plateau at high frequencies.For the treated CNT suspensions studied in this paper, a single-mode Maxwell model, however, could not describe the experimental LVE data which showed GЈ scaling roughly as 0.5 over three decades of frequency.To reconcile the difference between prediction by Eq. ͑25͒ and experimental LVE data, one might consider the concept of an "effective rotary diffusion coefficient."Most notably, Folgar and Tucker ͑1984͒ proposed the use of an empirical diffusion coefficient where D ˆr = C I ␥ ˙ϰ C to account for fiber-fiber interactions in a semidilute short fiber suspension ͑C is known as the Folgar-Tucker constant and was determined to be about 0.003-0.016͓͒Larson ͑1999͔͒.If it is assumed that D ˆr = C, this would give a slope of 1 on a log-log plot of GЈ versus frequency, inconsistent with the experimental results obtained for the treated CNT suspensions.It is believed that the physics of chemically treated CNTs is more complex and that the LVE rheology of treated CNT involves hydrodynamic interactions, Brownian motion, fiber-fiber collisions, as well as other possible effects, such as bending and repulsion, not considered by Folgar and Tucker ͑1984͒.
In this paper, we assume that the effective diffusion coefficient can be written as where is the base diffusion coefficient accounting for randomizing effects, such as Brownian motion and tube-tube interaction, and was determined to be 0.005 s −1 in fitting the orientation model to steady shear experimental data; D r2 is the diffusion coefficient corresponding to a weak network of CNT that is only present at small strain.͑␥͒ is a function which approaches 1 at a very small strain and equals to 0 for an applied strain ͑␥͒ larger than a critical strain ͑␥ c ͒ that destroys the network.
Replacing the diffusion coefficient D r term in Eq. ͑12͒ with the effective diffusion coefficient defined in Eq. ͑27͒ gives Phenomenologically, if it is assumed that D r2 scales as the square root of frequency ͑i.e., D r2 = C ͱ ͒, reasonable fitting to experimental LVE data can be obtained.The best fits for the 0.2% and 0.3% treated CNT suspensions are shown in Fig. 6.In these fittings, C = 0.05 s −0.5 ; whereas the concentration-dependent parameters N p and ␤ scaled with the CNT concentration.Although experimental GЈ and GЉ data were not available at frequencies lower than 0.1 s −1 ͑due to the small torque generated͒, the model predicted GЈ ϰ 2 at very low frequencies and this is in agreement with the Brownian rod models described, for example, in Larson ͑1999͒.The only difference between the current model and the Brownian rod model is that the latter model predicted a plateau at high frequencies which was not observed experimentally for the treated CNT suspensions.Mathematically, the current model reduced to the rod model by setting C = 0 as shown in Fig. 7.
Figures 7 and 8 show model predictions for GЈ, GЉ and ‫ء‬ by varying the values of C and ␤ respectively.It is clear from these figures that GЈ is very sensitive to the values of C and ␤ while GЉ is not.The results shown in Fig. 6 were obtained by combining Eqs.͑25͒ and ͑26͒ with the assumption that D r2 scaled with the square root of the frequency.The same results can be obtained by applying the same assumption and solving Eqs.͑12͒ and ͑13͒ with a linear closure relation and an effective diffusion coefficient given by Eq. ͑27͒.The consistency of the results confirms the approximations introduced in the derivation of Eqs.͑25͒ and ͑26͒ were appropriate.
Despite the success of the fitting, an empirical relation where D r2 = C ͱ was assumed and the exact physical meaning behind such a relation is unclear at this stage.It is, however, worth noting in the context of polymer kinetic theories that similar evolution of G គ Ј has been predicted by the Rouse model.Rouse model predicted that at low frequen- cies, GЈ ϰ 2 and at higher frequencies, GЈ ϰ 0.5 .The factor of 0.5 is a result of the spacing of relaxation times in the Rouse model ͓Bird et al. ͑1977͒; Larson ͑1999͔͒ and following a similar logic, it is probable that the experimental evolution of GЈ for chemically treated CNT suspensions was a result of multiple relaxation times.Experimental LVE data can be fitted using the Rouse approach by assuming CNTs as curved and elongated entities that can be divided into subsections; each represented by a pair of bead and spring.Such an approach, however, also implies that the elasticity of CNTs originated from the stretching of CNTs.At present, there is no consensus over the origin of mild elasticity for chemically treated CNT suspensions.Xu et al. ͑2005͒ proposed that the elasticity was due to stretching and bending of CNFs, whereas Hough et al. ͑2004͒ suggested that the elasticity originated from the interaction between CNTs instead of bending or stretching.The latter belief was also held by Dinsmore et al. ͑2006͒ in their treatment of colloidal dispersion, where they considered that the interactions between particles could be represented by conceptual springs.
In the current modeling, the treated CNTs were assumed to be essentially rigid, but we do not completely rule out the possibility that CNTs can bend in certain shear conditions.Duggal and Pasquali ͑2006͒ tagged some surfactant-stablized CNTs with fluorescent dye and were able to estimate the persistence length of CNT suspended in water.They computed the persistence length to be in the order of 7 -74 m ͑in the absence of any external field͒.Given the typical length of the treated CNTs in the current study was less than 1 m, we tend to believe that the treated CNTs were essentially rigid.For future work, it is of interest to study in greater detail the bending stiffness of the chemically treated CNT and identify the effect of bending ͑if any͒ on the rheological response of treated CNT suspensions.Intuitively, the elastic energy can be stored up due to the bending of CNT and could therefore potentially contribute toward the elasticity of the suspension.
Based on LVE experimental data alone, the exact origin for the mild elasticity of treated CNT suspensions remains inconclusive.Both explanations are possible and are not mutually exclusive.However, there is an evidence in step strain experiments showing the presence of a weakly interconnected CNT network.Therefore, we tend to believe that the weak network has contributed to the mild elasticity of treated CNT suspensions.

V. STEP STRAIN EXPERIMENTAL DATA
A series of step strain experiments were carried out in order to reveal more detailed relaxation behavior of the treated CNT suspensions and offer insights into the origin of elasticity.A finite step strain ͑␥ 0 ͒ was applied to the CNT suspensions and the process of stress relaxation was followed using the ARES strain-controlled rheometer.Figure 9͑a͒ shows the time evolution of the relaxation modulus ͑G͒, which is defined as G͑␥ , t͒ = / ␥ 0 , for suspensions with different concentrations of CNT.The stepper motor had a response time of about 0.1s ͑as indicated in the figure͒ and for the epoxy matrix, the stress dissipated almost instantaneously consistent with the fact that it behaved essentially as a simple Newtonian fluid in both steady shear and LVE experiments.Addition of CNTs prolonged the stress relaxation process and the CNT suspensions resembling the response of a viscoelastic fluid.The effect was progressive as the CNT concentration increased and this confirmed the earlier LVE experiments that the addition of CNTs effectively increased the elasticity of the system as a whole.
Different magnitudes of strain were applied to the 0.5% CNT suspension.Figure 9͑b͒ shows a strain dependence in terms of the final mode of stress dissipation.At small strains ͑1%, 5%, and 10%͒, the CNT suspension responded essentially as an entangled gel with a strain softening characteristic associated with a long relaxation time ͓see, for example, Mours and Winter ͑1996͒ and Winter ͑1999͔͒.At high strain, the CNT suspension dissipated in a dominantly viscous fluid manner.Similar strain-dependence phenomenon was reported by Guskey and Winter ͑1991͒ for a thermotropic liquid crystalline polymer in the nematic state.Intuitively, the strain-dependence relaxation process can be explained by the yielding of a network ͓Mewis and Meire ͑1984͒; Amari and Watanabe ͑1980͔͒.Depending on the strength of the network, if a large enough strain is applied, the network will be broken down and will finally dissipate as a fluid.The network for the 0.5% suspension is considered to be a relatively weak one and broke down at a strain level higher than 10%.The findings have two implications.First, it is highly probable that the mild elasticity as observed in LVE measurements is linked to the presence of a weak CNT network, and second, the effect of elasticity is negligible at high strain level, in line with the concept of a nonlinear damping parameter proposed by Wagner ͑1976͒.
In terms of the orientation model, the first term in Eq. ͑12͒ involving the strain rate tensor vanishes in the absence of flow after the application of step strain and the equation reduces to 12 = ␤D r a 12 .͑30͒ Similarly, Eq. ͑13͒ becomes If D r is assumed to be independent of a 12 , integration of Eq. ͑31͒ would give where C is an integration constant corresponding to the initial orientation.This gives 12 = C␤D r e −6D r t or G͑t͒ = 12 ␥ 0 = C␤ ␥ 0 D r e −6D r t .͑33͒ Equation ͑33͒ predicts an exponential decay in stress relaxation after the application of step strain.This is, however, inconsistent with the experimental results shown in Fig. 9 where relaxation modulus G͑t͒ was found to exhibit essentially a power-law decay over intermediate time scales ͑0.1 s Ͻ t Ͻ 10 s͒.Successful quantitative fitting to experimental data requires the D r in the orientation model to be a certain function of the orientation a 12 , as in the case of liquid crystal polymer modeling where the nematic potential depends on the molecular orientation distribution ͓see, for example, Doi and Edwards ͑1986͔͒.The orientation dependence of rotary diffusion coefficient could possibly be explained by the presence of repulsion force between CNTs after chemical treatment.This again supports the use of an effective diffusion coefficient term, which contains not only rotary Brownian diffusion but also other possible randomizing events such as repulsion.It is conjectured that CNTs were isotropically oriented prior to the step strain and subsequent step strain has led to the ͑partial͒ alignment of CNTs.In the absence of flow, the orientation of the partially oriented CNTs then began to randomize as a result of the Brownian motion and possible repulsion.Repulsion ͑if any͒ would be minimized when the CNTs were isotropically oriented and with their relative distance maximized; whereas the progressive alignment of CNTs would decrease the relative distance and therefore increase the repulsion between CNTs.This implies that the effective diffusion coefficient would start off at a high value after the initial step strain oriented some of the CNTs and its magnitude would decrease as the orientation of CNT became randomized again.It is clear from the analysis that quantitative fitting to step strain experimental data is possible, but it also requires certain assumption about D r as a function of a 12 , further increasing the complexity of the current model.For this reason, the quantitative fitting to step strain experimental data is not included in this paper, of which the main focus is to develop a self-consistent model for describing both steady shear and LVE data.

VI. A GENERALIZED ORIENTATION MODEL FOR DESCRIBING BOTH STEADY SHEAR AND LINEAR VISCOELASTICITY DATA
In Sec.III, a FP based model was successfully developed for steady shear flow using a rotary diffusion coefficient D r1 as a primary variable parameter.In Sec.IV, again a FP based equation can be derived to reproduce the LVE response and in this case using a differently formulated diffusion coefficient D r2 .Clearly, it is desirable to have a constitutive equation that is self-consistent for general deformation, including steady shear and LVE.The situation has similarities to the viscometric response of molten polymers, where Maxwell modeling successfully described LVE and an additional strain softening parameter was introduced to describe the full nonlinear response ͓Wagner ͑1976͔͒.
For the case of the FP formulation, it is possible to reconcile both LVE and nonlinear steady shear response by formulating the constitutive equation in the following way: Constitutive equation: where In general, the second term in the constitutive equation becomes negligible compared with the first term in steady shear and large strain deformation ͓Eqs.͑15͒ and ͑16͔͒ and the effect of randomizing events is mainly accounted for by the FP equation.In both the constitutive equation and FP equation, the rotary diffusion coefficient is replaced by an effective diffusion coefficient D ˆr = ͑D r1 + ͑␥͒D r2 ͒.The effective diffusion coefficient consists of two elements: ͑1͒ a constant diffusion coefficient D r1 representing Brownian motion and tube-tube hydrodynamic interactions.Its value can be identified from the steady shear data fitting and a constant value was found to be sufficient in describing steady shear responses; ͑2͒ an empirical diffusion coefficient D r2 that corresponds to a weakly interconnected CNT network and it is active only at small strains.Phenomenologically, D r2 was found to be proportional to the square root of frequency in smallamplitude oscillatory flows.When the strain applied exceeds a certain critical strain ͑␥ c ͒, the weak network is destroyed and only D r1 remains.

VII. CONCLUSIONS
Experimental results for shear thinning and viscoelastic data have been presented for treated CNT suspensions.A FP orientation diffusion model has been developed and has been shown to predict the base shear-thinning response for certain chemically treated CNT suspensions.Untreated CNT suspensions in epoxy resin have been reported to exhibit an optical microstructure ͓Rahatekar et al. ͑2006͒; Ma et al. ͑2008a͔͒; however, treated CNT suspensions reported in this paper showed no discernable optical microstructure, indicating that the suspensions were well dispersed at least on a micronscale.The modeling of treated CNTs involves solving the FP diffusion equation in a closed form and fitting the experimental shear-thinning data to two adjustable parameters N p and D r1 .The quality of the fit was considered reasonable although a single value of rotary diffusion coefficient D r1 was insufficient to fit the whole range of data exactly.The success of the model suggests that the CNTs are behaving essentially as high aspect ratio rigid rods with a very low rotary diffusion coefficient and that shear is able to align the rods in the direction of flow.The findings have relevance to the processing of CNTs where in some cases a preferential orientation can have a strategic advantage and there are other situations where an isotropic distribution of CNTs is desired.
In terms of LVE properties, CNT suspensions exhibited an elastic behavior that is different from ordinary non-Brownian fiber or Brownian rod suspensions.Using a rotary diffusion term that does not depend on the applied frequency was found to be inappropriate in describing the evolution of storage modulus ͑GЈ͒ and loss modulus ͑GЉ͒ as a function of frequency.The experimentally observed elasticity was considered to be relatively mild and suppressed when subject to large deformations, according to step strain experiments.It is conjectured that the mild elasticity is due to the presence of a weakly interconnected CNT network.Experimental GЈ and GЉ data were fitted using an effective diffusion coefficient and an empirical relation was identified.Finally, a unified FP based orientation model was formulated to describe both steady shear and viscoelastic responses of the treated CNT suspensions.which leads to a linear system of N ϫ ͑M −2͒ algebraic equations involving N ϫ M nodal values i,j n .This system of equations must be completed with the equations for the upper and lower nodes ͑j = 1 and j = M͒.At these nodes, a numerical difficulty appears because all the upper nodes ͑i, j = M͒ correspond to the same point on the unit sphere ͑the upper pole͒ and the same for the lower ones.Due to the symmetry of the distribution function, the same unknown value is assigned to all these nodes, that is, The N ϫ M discrete equations involve N ϫ M + 1 unknowns, N ϫ M nodal values i,j n , and the introduced value.An additional equation must be introduced to define the solution.For this purpose, the discrete form of the normality condition ͓Eq.͑5͔͒ can be introduced, which leads to where A ij represents an area on the unit sphere that is related to the node ͑i , j͒, A ij = h sin͑ j ͒h .͑A13͒ The linear system can then be solved to obtain the nodal values of the orientation distribution function i,j n at time t = n⌬t.The discretisation of the Fokker-Planck equation leads to a system of algebraic linear equations and each equation represents the evolution of a fraction of tubes oriented in a certain direction.Each of these linear equations corresponds to a particular orientation on the unit sphere and represents the balance between the tubes coming from and leaving to its "neighboring" orientations, as a result of shear flow and other randomizing events such as Brownian motion and tube-tube hydrodynamic interaction.It was verified that the numerical solutions reported in this manuscript did not depend on the numbers N and M, further confirming the solutions were good approximation to the actual solutions of the kinetic equation Eq. ͑8͒.
The time step must be carefully chosen to ensure the stability conditions related to the evolution problem and the convective character.If the time step becomes too small due to large Peclet numbers defined at high shear rates, some upwinding stabilization could be applied.͑The Peclet number is defined as the ratio between the maximum fiber rotation velocity and the diffusion coefficient.͒However, it is well known that the upwinding introduces certain amount of nonphysical over diffusion.To avoid these nonphysical effects, simulations were performed without any stabilization, although more computation time is required.
FIG. 4. ͑a͒Orientation model fitting for epoxy, 0.05%, 0.2%, 0.33% and 0.5% chemically treated CNTs suspended in epoxy where N p and D r are model fitting parameters.͑b͒ Sensitivity of the orientation model to the fitting parameter D r for the 0.33% treated CNT suspension and N p =7. ͑c͒ The scalar order parameter ͑S͒ calculated by fitting apparent shear viscosity data to the orientation model for a rotary diffusion coefficient of 0.005 s −1 .͓Note: by setting N p = 0, the fiber contribution term in Eq. ͑16͒ becomes zero and therefore, a = where is the suspending medium viscosity.͔

͑19͒
FIG. 5. ͑a͒Loss modulus ͑GЉ͒ and ͑b͒ storage modulus ͑GЈ͒ as a function of frequency for different concentrations of treated CNT suspensions ͑with strain= 1%͒.͑c͒ LVE data, which include the storage modulus ͑GЈ͒, the loss modulus ͑GЉ͒, and the complex viscosity ͑ ‫ء‬ ͒ as a function of frequency, for the 0.5% treated CNT suspension.

FIG. 6 .
FIG. 6.Model fits to experimental LVE data of ͑a͒ 0.2% and ͑b͒ 0.5% treated CNTs in epoxy.N p , ␤, and C are fitting parameters and the concentration-related parameters N p and ␤ scale with the concentration.
and nanodevices ͓Ajayan et al. ͑1994͒; Calvert ͑1999͒; de Heer et al. ͑1995͒; Rinzler et al. ͑1995͒; Saito ͑1997͒; Collins et al. ͑2000͒; Kong et al. ͑2000͒; Tans et al. ͑1998͔͒.Most of these applications involve suspending the CNTs in a matrix and thereby allows for further processing such as fiber spinning, film casting, extrusion, and inkjet printing ͓Vigolo et al. ͑2000͒; Safadi et al. ͑2002͒; Ericson et al. ͑2004͒; Kozlov et al. ͑2005͒; Kordás et al. ͑2006͔͒.Both in terms of processing and product, a uniformly dispersed suspension of CNT is generally desirable ͓Shaffer et al. ͑1998͒; Sandler et al. ͑1999͔͒; however, CNTs have low solubility in most common solvents and therefore obtaining a homogenous dispersion of CNT can be laborious and difficult ͓Calvert ͑1999͒; Dresselhaus and Dai ͑2004͔͒.For instance, Chaubal et al. ͑1997͒; Chinesta et al. ͑2003͒; Lozinski and Chauviere ͑2003͒; Chauviere and Lozinski ͑2004͒; Keunings ͑2004͒; Ammar and Chinesta ͑2005͒; Ammar et al. ͑2006a, 2006b͔͒.In the subsequent LVE modeling, given that the orientation distribution of CNT during small-amplitude oscillations was expected to remain close to an initial isotropic distribution, a linear closure relation was used.