Amine Ammar
email: amine.ammar@ujf-grenoble.fr

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Elías Cueto
email: ecueto@unizar.es

Coupling finite elements and proper generalized decompositions

Keywords: Coupling finite elements and proper generalized decompositions Separated representations, Finite sums decomposition, Finite Elements, Enriched Finite Elements, Proper Generalized Decomposition

Introduction

Numerous models encountered in science and engineering remain nowadays, despite the impressive recent progresses attained in computational simulation techniques, intractable when the usual and well experienced discretization techniques are applied for their numerical simulation. Thus, different challenging issues are waiting for the proposal of new alternative advanced simulation techniques, the brut force approach being no more a valuable alternative.

A first challenging issue concerns the treatment of highly multidimensional models arising from quantum mechanics or kinetic theory descriptions of solids and complex fluids, including micro and nano-structured materials [START_REF] Chinesta | The nanometric and micrometric scales of the structure and mechanics of materials revisited: An introduction to the challenges of fully deterministic numerical descriptions[END_REF]. Other multidimensional models are encountered in biology, genetics and chemistry where the so called chemical master equation plays a key role [START_REF] Sreenath | Modelling the dynamics of signalling pathways[END_REF]. Curse of dimensionality also arises in stochastic models involving numerous random variables [START_REF] Nouy | Recent developments in spectral stochastic methods for the solution of stochastic partial differential equations[END_REF]. The main challenge in the treatment of this kind of models is related to their multidimensionality because when one applies standard mesh based discretization techniques the number of degrees of freedom involved scales exponentially with the dimension of the space concerned. Thus, in high dimensions, usual meshes cannot be defined at all, needing the proposal of new advanced strategies able to circumvent the terrible curse of dimensionality. Until now, the treatment of this kind of models was restricted to the ones defined in moderate multidimensional spaces where for example the sparse grid based methods works [START_REF] Bungartz | Sparse grids[END_REF].

We proposed in some of our former works [START_REF] Ammar | A new family of solvers for some clases of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF] [2] a separated representation strategy able to circumvent the curse of dimensionality that those models involve. Thus, given a model involving the unknown field u(

x 1 , x 2 , • • • , x D ), with (x 1 , x 2 , • • • , x D ) ∈ Ω 1 × Ω 2 × • • • × Ω D , where Ω i ⊂ R di , ∀i, 1 ≤ i ≤ D, the solution is searched as u(x 1 , x 2 , • • • , x D ) ≈ i=N i=1 α i F i 1 (x 1 ) • • • F i D (x D ) (1)
where the dimension of the model is DIM = j=D j=1 d j . Obviously, if a mesh is defined in each domain Ω i consisting of N i n nodes, the whole mesh of Ω = Ω 1 × • • • × Ω D involves i=D i=1 N i n nodes. We can notice that the number of degrees of freedom increases exponentially with the dimension of the space. On the contrary, by assuming the separated representation (1) the number of degrees of freedom reduces to: N × i=D i=1 N i n proving that the complexity scales linearly with the dimension of the space.

Another issue concerns the solution of transient multiscale models (usually strongly non linear and coupled, and always of high size). These models arise in computational mechanics (involving a large variety of constitutive behaviors, couplings etc.). In this context, the use of standard incremental discretization techniques becomes inefficient from the computational time viewpoint. Moreover, in general, these models involve different characteristic times differing of several decades. Again, alternative advanced computational techniques are needed. Pierre Ladeveze proposed several years ago a powerful technique for addressing this kind of challenging models that he called the LATIN method [START_REF] Ladeveze | Non linear computational structural mechanics[END_REF].

The LATIN method integrates many ingredients leading to a robust, powerful, efficient and accurate discretization technique especially well adapted for treating transient multi-scale non-linear models. The two most outstanding ingredients are (i) the decoupling between a linear-global problem and a nonlinear-local one, both defined in the whole space-time domain; and (ii) a spacetime separated representation of the model variables in order to accelerate the solution of the linear-global problem. The former separated representation was called by Ladeveze in the 80's "radial approximation", and in our knowledge it was the first time that separated representations were applied in computational mechanics. The interested reader can refers to [START_REF] Ladeveze | The LATIN multiscale computational method and the proper orthogonal decomposition[END_REF] and the references therein for a recent overview and the state of the art of the LATIN approach in multi-scale modeling.

Thus, the radial approximation proposed by Ladeveze can be viewed as a particular case of Eq. ( 1) with only two coordinates, i.e. u(x 1 , x 2 ), x 1 representing the physical space, i.e. x 1 ≡ x ∈ Ω ⊂ R d (d = 2 or d = 3); and x 2 the time, i.e. x 2 ≡ t ∈ I ⊂ R + . However, because the radial approximation was proposed before the separated representation given by Eq. ( 1) (originally proposed in [START_REF] Ammar | A new family of solvers for some clases of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF] [2]) we could also state that Eq. ( 1) is no more than a generalization of the so called radial approximation for models involving many coordinates.

Separated representations were applied before its use in computational mechanics in other branches of the computational physics. A remarkable example concerns quantum chemistry that focuses in the modeling of atomic structure and molecular bonds. In that case, one must solve the Schrödinger equation, or some more advanced alternatives, as the fully relativistic Dirac equation. The solution of these equations is of capital interest to describe the structure and mechanics of materials at the finest scale, descriptions that can be directly applied in the framework of nanosciences and nanotechnologies, but that also allow deriving coarse grained models by applying appropriate multiscale bridges. Thus, quantum chemistry leads to "ab initio" simulations that at its turn leads to molecular dynamics simulations. Coarse graining results first in Brownian dynamics simulations and then on statistical mechanics descriptions. The coarsest scale concerns obviously the macroscopic modeling. This multiscale modeling needs accurate models and solvers at each scale as well as appropriate bridges for moving from one scale to the neighbor ones.

However, the solution of the quantum scale models is quite delicate because the wavefunction that describe for example the electronic distribution, is defined in a highly multidimensional space whose dimension scales with the number of particles that the system involves. Thus, the only possibility of representing a function in such a space was the use of separated representations that constitutes the basis of the Hartree-Fock and the post-Hartree-Fock approaches [START_REF] Cancès | Computational Quantum Chemistry: a primer[END_REF]. Despite the nice dream that this multiscale modeling (ranging from the quantum to the macroscopic scale) opened, this procedure is in fact inapplicable at present, and probably for many many years, because the starting point fails: the Schrödinger or the Dirac equations remain intractable due to the Pauli's exclusion principle constraint. Simplified models are then mandatory (see [START_REF] Cancès | Computational Quantum Chemistry: a primer[END_REF] for an excellent overview).

Numerous topics concerning the use of separated representation were and are being addressed in our researches. Some recent developments concern:

• From the mathematical point of view some results concerning the convergence were obtained, but nothing is known about the rate of convergence [START_REF] Ammar | A Kronecker product based iterative solvers for linear systems, Numerical Linear Algebra with Applications[END_REF]. Its connection with the POD (Proper Orthogonal Decomposition) was also analyzed in [START_REF] Ladeveze | The LATIN multiscale computational method and the proper orthogonal decomposition[END_REF]. In fact, introducing the separated representation into the weak formulation of the model the different functions involved in the separated representation can be computed. However, these functions that are optimal from the point of view of the strategy considered, are no more orthogonal. For this reason, we called the resulting approximation Proper Generalized Decomposition (PGD).

• From the point of view of the applications, there is a huge catalogue [START_REF] Nouy | Multiscale computational strategy with time and space homogenization: A radial-type approximation technique for solving microproblems[END_REF] [1] [2] [3] [8] among many others.

• From the point of view of the computing time savings, these are often of some orders of magnitude and in some cases this technique allowed solving models never until now solved, implying hundreds of dimensions, in the context of deterministic or stochastic models [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: Transient simulation using space-time separated representation[END_REF].

• We proved recently that in some non-linear-local models arising from the natural decoupling of a non-linear-global model (the LATIN strategy being the most natural choice) can be globalized by applying a separated representation. Even for strongly non-linear behaviours (e.g. thermal models involving strongly non-linear curing kinetics usually encountered in polymer or composites forming processes) the CPU time savings are impressive.

• Some proposals concerning the adaptivity and the error estimation have been done. However, further developments are needed to improve their efficiency.

• The stabilization of models involving advective terms is another issue that needs further developments. Another open issue concerns the stabilizations of mixed formulations for verifying the stability (LBB) conditions.

• Separated representations of the 3D physical space can be performed, and easily coupled with standard or enriched finite element descriptions for capturing evolving discontinuities, singularities or to perform subdomains coupling.

This paper focuses in one of the just referred topics, the one concerning the separated representation of variables defined in usual physical spaces (2D or 3D). In [START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF] we applied such one approach to perform high resolution homogenization where the representative volume (a 3D cube in that case) was described using thousand nodes along each of the coordinate directions. In that case, instead of using a mesh involving 10 9 degrees of freedom, we considered a separated representations that only involved few thousands of degrees of freedom. At present other groups are considering this separated representation to perform direct numerical simulation in flows involving turbulence, identification of permeability in porous media and homogenization in heterogenous media.

Recently, we extended the use of separated representations to domains others that cubes or hyper-cubes where the issue of non-homogeneous boundary conditions were also addressed [START_REF] Gonzalez | Recent advances in the use of separated representations[END_REF]. However, models can also involve interfaces where the field or some of its derivatives could be discontinuous. In this case the use of separated representation seems to be delicate. We prove in this work that a coupling between the separated representation for describing the smooth solution behavior and a finite element (or an enriched finite element) description in some patches of the domain where the solution is expected to exhibit localization or some kind of discontinuity, seems to be an appealing choice to address these models.

Other issue concerns the coupling of different subdomains, each one involving a different separated representation of the solution. Thus, the required continuity conditions at those interfaces must be enforced. The present work analyzes this difficulty and propose an efficient strategy for coupling subdomains also based on the use of appropriate finite element patches. A direct consequence of this last procedure is the coupling of two subdomains, one of them involving a separated representation of the unknown fields, whereas the other one considers a finite element description. Thus, one could easily couple the finite element technology, well experienced and very powerful, with a reduced modeling performed by using separated representations in the regions of the domain where the solution is expected evolving smoothly (in many cases these regions represent the largest part of the domain in which the model is defined, as in the case of many models in fracture mechanics where the evolving crack is very localized and only affects a small region of the material domain; models in fluid-structure interaction, where that interaction only occurs in a very localized region of a larger domain, ...)

Next section revisits the algorithm related to the PGD of the unknown field. In section 3 this strategy is generalized by considering a finite element based enrichment of the separated representation for capturing localized behavior. Finally, section 4 describes some potential applications of this enriched formulation, in particular a model that exhibits localization and another one involving subdomains coupling.

The Proper Generalized Decomposition revisited

For the sake of simplicity we are considering a simple 2D model in wich the unknown field depends on the coordinates x and y, defined in the domains Ω x and Ω y respectively. Now, the solution is searched in the form:

u(x, y) ≈ i=N i=1 α i F i (x) • G i (y) (2)
In order to build-up this separated representation, an iterative algorithm is proposed, solving at each iteration a projection and an enrichment problem.

If we assume that the first n functions in the sum has been already computed

u(x, y) ≈ i=n i=1 α i F i (x) • G i (y) (3)
the projection stage consists of finding the n best alpha-coefficients. Then, the enrichment stage lies in finding the best functions R(x) and S(y) such that the updated representation given by

u(x, y) ≈ i=n i=1 α i F i (x) • G i (y) + R(x) • S(y) (4)
satisfies the weak formulation for particular test functions. We are assuming a generic equation

K (u(x, y)) + L (u(x, y)) = 0 ( 5
)
where K and L are two differential operators. For the sake of simplicity we are assuming that the first one only involves derivatives with respect to the x-coordinate, the second one involving the derivatives with respect to the other coordinate.

The computer implementation of these steps needs for a discrete representation of all the functions: F i (x), G i (y), R(x) and S(y) using a finite element interpolation in their associated domains Ω x and Ω y . For this purpose we introduce the vectors M and N containing the shape functions associated with the meshes of Ω x and Ω y respectively. Finally F i , G i , R and S represent the nodal description of the associated functions. We define the following matrices related to the weak formulation of equation ( 5):

         M = Ωx MM T dΩ x N = Ωy NN T dΩ y K = Ωx MKM T dΩ x L = Ωy NLN T dΩ y (6)
These integrals take into account the specific character of each operator. For example integration by parts is used in second order operators, upwinding for stabilizing advective terms, etc.

Projection stage

In this stage we consider the test functions given by

u * (x, y) = n i=1 α i * F i (x) • G i (y) (7)
which writes in the discrete form:

u * (x, y) = n i=1 α i * F iT MG iT N (8)
Using the discrete form of function u(x, y):

u(x, y) = n i=1 α i M T F i N T G i (9)
the weak formulation of Eq. ( 5) leads after numerical integration to:

n i=1 n j=1 α i * H ij α j = 0 (10)
where

H ij = F iT KF j • G iT NG j + F iT MF j • G iT LG j (11)
This system must be solved taking into account boundary and/or normality conditions. For example, we could consider a first term in the sums decomposition verifying the boundary conditions, and then the other ones must vanish on the domain boundary (this condition is enforced in the enrichment stage). In this case we must enforce α 1 = 1.

Enrichment stage

In this stage we consider the unknown field given by

u(x, y) = n i=1 α i F i (x) • G i (y) + R(x) • S(y) (12)
and the test functions by

u * (x, y) = R * (x) • S(y) + R(x) • S * (y) (13)
whose discrete form writes

u * (x, y) = R * T S * T M(S T N) (R T M)N (14)
The weak form of Eq. ( 5) becomes, in the framework of a fixed point strategy, in:

• For a given S(y) (that implies u * (x, y) = R * (x) • S(y))

KS T NS + MS T LS R = = - n i=1 α i KF i S T NG i + MF i S T LG i (15)
• For a given R(x) (that implies u * (x, y) = R(x) • S * (y))

R T KRN + R T MRL S = = - n i=1 α i R T KF i NG i + R T MF i LG i (16)
It must be noticed that the number of the degrees of freedom involved in such non linear solution is the sum of the degrees of freedom involved in each direction (instead the product of them that mesh based discretizations strategies would have involved).

Finally, after convergence, the new approximation functions F n+1 and G n+1 are obtained by normalizing R and S:

F n+1 = R(R T MR) -1 2 G n+1 = S(S T NS) -1 2 (17)
The projection and enrichment stages must continue until reaching convergence.

Finite elements based separated representations enrichment

In this section we are modifying the technique just presented in order to capture localized behaviors such as high gradients or weak and strong discontinuities. This kind of behaviors induces too many functional products within a standard separated representation. However, it is well known that they can be easily taken into account within the finite element framework, by using adaptive local mesh refinement in the first case or any of the techniques able to represent weak and strong discontinuities (meshes compatible with the interfaces or enriched approximations when the meshes are not compatible with the interfaces).

Thus, the mail goal is how combining separated representations, able to approximate the solution in the most part of the domain where the model is defined, reducing significantly the computational cost, and a localized finite element description, that even if it is expensive in nature, it only applies in a small region of the whole domain where the solution is expected to exhibit localized behavior. To this end, we propose a multi-scale approach that somewhat resembles the s-version of the finite element method by J. Fish [START_REF] Fish | The s-version of the finite element method for laminated composites[END_REF] or the multi-scale FEM proposed by Rank [START_REF] Rank | A multiscale finite element method[END_REF], but in this case in a global (Ritz) and separated basis function setting.

We come back to the model described in the previous section (Eq. ( 5)) defined in Ω = Ω x × Ω y and we imagine that the solution u(x, y) is smooth enough everywhere except in a small region Ω l ⊂ Ω. Now, we could imagine an approximation combining both, the separated representation and a finite element approximation, where the last one only applies in Ω l . We assume a mesh on Ω l composed of N l nodes. We also assume that the finite element shape functions related to the finite element approximation cannot be expressed from the tensor product of the one dimensional bases employed to build-up the separated representation according to the description given in the previous section.

Now, the approximation in the whole domain can be written as:

u(x, y) = u SR (x, y) + u F E (x, y) (18)
where u SR (x, y) is defined in the whole domain Ω = Ω x × Ω y whereas the finite element enrichment u F E (x, y) is only defined within Ω l vanishing in Ω -Ω l . In order to ensure the continuity of the resulting approximation we must enforce the nullity of the enrichment u F E (x, y) on the boundary of Ω l , ∂Ω l . The resulting approximation writes:

u(x, y) ≈ i=N i=1 α i F i (x) • G i (y) + j=N l j=1 Ψ j (x, y)u j (19)
where Ψ j (x, y) are the standard finite element shape functions and u j the associated weights. Because the contribution of u SR (x, y) in Ω l , u j do not correspond to the values of the unknown field u(x, y) at the nodal positions (x j , y j ).

Remark.-We are assuming that the enrichment is performed by using standard finite elements, but in fact any kind of compact support approximation could be used. In order to represent interfaces involving weak discontinuities one could proceed within the standard finite element method by ensuring that the interface coincides with the elements edges. If the interface passes across the elements an enriched version of the finite element method should be considered (e.g. the extended finite element [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF]). Now, to build-up the just defined approximation [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF] we proceed as in the previous section by alternating a projection and an enrichment step until reaching convergence. The degrees of freedom related to the finite element contribution u j will be computed at the projection step because at this step we consider all the approximation functions known and we look for the associated weights. At the enrichment step a new functional couple of the separated representation will be calculated by assuming known all the previously computed functional couples as well as the just updated finite element enrichment. We are summarizing both steps.

Projection stage

At this stage we consider the test functions given by

u * (x, y) = n i=1 α i * F i (x) • G i (y) + N l j=1 Ψ j (x, y)u * j (20)
which writes in the discrete form:

u * (x, y) = n i=1 α i * F iT MG iT N + U * T Ψ ( 21
)
where Ψ is the vector that contains the finite element shape functions related to the set of nodes considered within the enrichment region Ω l , and U the vector that contains the nodal degrees of freedom.

Using the discrete form of function u(x, y):

u(x, y) = n i=1 α i M T F i N T G i + Ψ T U (22)
the weak formulation of Eq. ( 5) leads after numerical integration to the linear system

(α * T U * T ) H B C Q α U = 0 (23)
In Eq. ( 23) matrix H is coming from the contribution of u SR * and u SR and therefore it results the same that the one concerned in Eq. [START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF].

On the other hand Q concerns the contribution coming from u F E * and u F E that represents a classical finite element discretization in Ω l . The remaining terms B and C represent the crossing contributions.

In order to give an explicit form of the different matrices involved in Eq. ( 23), we introduce the following notation

         M SR = Ωx MM T dΩ x N SR = Ωy NN T dΩ y K SR = Ωx MKM T dΩ x L SR = Ωy NLN T dΩ y (24) K F E = Ω l ΨKΨ T dΩ L F E = Ω l ΨLΨ T dΩ (25) K SR-F E j = Ω l MN T KΨ T j dΩ L SR-F E j = Ω l MN T LΨ T j dΩ (26)
and

K F E-SR j = Ω l Ψ j KMN T dΩ L F E-SR j = Ω l Ψ j MLN T dΩ (27)
where, as mentioned in the previous section, these integrals take into account the specific character of each operator, as for example integration by parts or upwinding for stabilizing advective terms. These integrals, when they involves finite elements are performed using standard Gauss quadrature schemes. In the previous expressions Ψ j denotes the shape function related to node j, that is, the j-component of vector Ψ.

Using these notations, the different matrices in Eq. ( 23) write:

H ij = F iT K SR F j • G iT N SR G j + F iT M SR F j • G iT L SR G j (28)
Q ij = K F E + L F E ij ( 29
)
B ij = F iT K SR-F E j + L SR-F E j G i (30) C ij = F jT K F E-SR i + L F E-SR i G j (31)
We can notice from the previous expressions that system (23) is nonsymmetric.

Eq. ( 23) must be solved again taking into account boundary and/or normality conditions.

Enrichment stage

At this stage we consider the unknown field given by

u(x, y) = n i=1 α i F i (x) • G i (y) + N l j=1 Ψ j (x, y) • u j + R(x) • S(y) (32)
where the first two sums in the right member are assumed known (computed at the projection stage).

The test function reads:

u * (x, y) = R * (x) • S(y) + R(x) • S * (y) ( 33
)
whose discrete form writes

u * (x, y) = R * T S * T M(S T N) (R T M)N (34)
The weak form of Eq. ( 5) becomes, in the framework of a fixed point strategy, in:

• For a given S(y) (that implies u * (x, y) = R * (x) • S(y))

K SR S T N SR S + M SR S T L SR S R = = - n i=1 α i K SR F i S T N SR G i + M SR F i S T L SR G i -H R-F E (35)
where the finite element enrichment contribution H R-F E writes:

H R-F E = j=N l j=1 K SR-F E j + L SR-F E j S u j (36) • For a given R(x) (that implies u * (x, y) = R(x) • S * (y)) R T K SR RN SR + R T M SR RL SR S = = - n i=1 α i R T K SR F i N SR G i + R T M SR F i L SR G i -H S-F E (37)
where the finite element enrichment contribution H S-F E writes:

H S-F E = j=N l j=1 K SR-F E j + L SR-F E j T R u j ( 38
)
It must be noticed that the number of the degrees of freedom involved in such non linear solution is the sum of the degrees of freedom involved in each one-dimensional discretization, instead the product of them that results in mesh based discretization strategies.

Finally, after convergence, the new approximation functions F n+1 and G n+1 are obtained by normalizing R and S:

F n+1 = R(R T MR) -1 2 G n+1 = S(S T NS) -1 2 (39)
The projection and enrichment stages must continue until reaching convergence. Two simple possibilities exist: (i) evaluating the norm of R(x) and S(y) and use these norms as stopping criterion (this criterion does not guarantee the convergence), or (ii) evaluate the norm of the strong formulation residual that is more expensive but it guarantees the convergence. In what follows we are using the second stopping criterion.

A first numerical example involving localization

To illustrate the capabilities of the just described strategy we are solving a simple Poisson equation

-∆u = s(x, y), (x, y) ∈ Ω = - π 2 , π 2 × - π 2 , π 2 (40)
with homogeneous boundary conditions, i.e.

u(x ∈ ∂Ω) = 0 (41)
The source term s(x, y) considered in the numerical experiment described in this section was:

s(x, y) = 1 + -4f (s 1 + s 2 s 3 ) s 4 (42)
with

       s 1 = -1 + 2d 2 f + f x 2 + f y 2 + 2df (x + y) s 2 = e 4df (x+y) s 3 = (-1 + 2d 2 f -2df (x + y) + f (x 2 + y 2 )) s 4 = e f (2d 2 +x 2 +y 2 +2d(x+y)) (43)
where coefficient d controls the location of the localized behavior whereas coefficient f controls the support of such localized behavior. In the numerical experiment that follows we considered f = 40 and d = 1.4. The above source term admits a separated representation involving 5 functional products.

The particular form of this source term leads to a localized behavior of the solution in the neighborhood of the upper-right and lower-left corners. This behavior allows to analyze the improvement that an eventual finite element based enrichment of the separated representation could bring.

First, the problem is solved using the standard finite element method. Then, it is solved again using a fully separated representation. Obviously, the existence of a localized behavior requires a fine enough resolution in both directions (one-dimensional meshes fine enough to capture the localized behavior) and also numerous functional products for approximating accurately the solution behavior.

Finally, the use of a separated representation combined with a finite element enrichment in the regions where the solution exhibits high gradients allows to define a technique sharing the benefits of both the finite elements and the separated representations.

Figure 1 depicts the mesh considered for solving problem (40)-(42) by using the finite element method. The associated solution is depicted in figure 2. Now, we are solving the same problem by using a coarse fully separated representation. The solution was computed by applying the technique described in section 2. We consider a very coarse separated representation, involving 7 nodes on each direction. As expected, the computed solution depicted in Fig. 3 does not account for the presence of the localized behavior. Figs. 4 depicts the one-dimensional functions associated with the coordinates x, F i (x) (on the left) and y, G i (y) (on the right), involved in the separated representation.

F i (x) (left) and G i (y) (right).
In this case, the regions in which the solution exhibits localization remain invisible to the separated approximation. Obviously, finer approximations of the one-dimensional approximation involved in the separated representation could finally capture the solution localized behavior. Instead of proceeding by refining the separated representation, our goal is capturing the localized behavior by superposing a finite element approximation in the regions exhibiting localization.

Figure 5 depicts the finite element meshes considered in the regions where the solution localization is expected.

Figure 5: Finite element meshes on the regions that exhibit localization. Now, we are using the finite element based enrichment of the separated representation, that is the mixed approximation that involves the separated form plus the contribution of the finite element approximation that is superposed to the former one in the regions depicted in Fig. 5. In this case the computed solution is shown in Fig. 6.

Obviously, as the one-dimensional meshes considered for defining the interpolations of functions F i (x) and G i (y) are too coarse to detect the presence of the localized behavior, the contribution of the separated representation to the whole solution depicted in Fig. 6 coincides to the one computed previously by using a fully separated representation (Fig. 3) and in consequence functions F i (x) and G i (y) involved in the enriched separated representation also coincide to the ones depicted in Fig. 4. Obviously, finer separated representations could represent a part of the localized behavior.

Obviously, the finite element enrichment captures the localized behavior as Fig. 7 depicts. This figure only represents the finite element contribution. We can notice that by construction, the finite element approximation vanishes on the boundary of the finite element meshes.

The coarse separated representation without the finite element representation involves a L 2 error of 0.23 (with respect to the finite element solution previously depicted). Obviously, the introduction of more terms in the finite sums decomposition cannot improve significantly this error because the onedimensional meshes cannot capture the localized behaviors as we commented previously. However, as soon as the finite element enrichment applies the error is drastically reduced to 0.09. For further error reductions finer finite element meshes should be considered. The use of finer separated representations could be also a strategy, but the effect of finite element remeshing is more efficient, because finite element interpolations are more adapted to capture localized behaviors.

A strategy for coupling subdomains

The enrichment procedure presented in this paper opens new possibilities for extending the applicability of separated representations. One of these potential applications lies in the coupling of subdomains, each one involving a separated approximation of the model variable. To illustrate the coupling procedure just described we are considering the solution of the model

-∆u = 2(y -y 2 ) + 2(1 -x 2 ), (x, y) ∈ Ω = (-1, 1) × (0, 1) (44)
with homogeneous boundary conditions, i.e.

u(x ∈ ∂Ω) = 0 (45)

The exact solution writes

u ex (x, y) = (1 -x 2 )(y -y 2 ) (46)
The domain Ω = (-1, 1) × (0, 1) can be decomposed in two subdomains Ω 1 = (-1, 0) × (0, 1) and Ω 2 = (0.1) × (0, 1) such that Ω = Ω 1 ∪ Ω 2 and Γ = Ω 1 ∩ Ω 2 (being Γ the interface between both subdomains). It is easy to understand that the solution of such one model is quite delicate because we should build-up a separated approximation of the solution in both subdomains assuring the appropriate transfer conditions across the interface Γ. Obviously, in the scenario just described the best alternative is the solution of the problem directly in the whole domain Ω because as it is rectangular, its decomposition in two subdomains is not needed. However, other kind of domains, e.g. L or T shaped domains, need appropriate subdomain decompositions.

An alternative lies in the construction of a separated representation in both subdomains assuming that both separated representations vanish on the interface Γ. Now, a finite element based enrichment is introduced in a strip Ω l of width 2H covering the interface, Ω l = (-H, +H) × (0, 1). As described in section 3, the nodal degrees of freedom related to nodes located on ∂Ω l are enforced to vanish, to ensure the compatibility (continuity) of the approximations across the enrichment domain boundary as well as for accounting the boundary conditions prescribed on ∂Ω. Thus, the finite element model defined in Ω l ensures the compatibility of the solutions computed in both subdomains in a finite element sense.

To perform the coupling of the model described above, we consider a strip covering the interface, Ω l = (-H, +H) × (0, 1), with H = 0.25, that is meshed using linear triangular finite elements. Fig. 8 shows the domain decomposition as well as the mesh of Ω l .

The exact solution is depicted in Fig. 9. Now, the solution in the whole domain is computed by applying the strategy just described. Figure 10 depicts the separated representation solution in Ω 1 (left) and Ω 2 (right) as well as the finite element enrichment in Ω l (center). We can see that the separated representation contribution in both domains vanish by construction on ∂Ω 1 and ∂Ω 2 . Thus, the separated representation contributions do not introduce any continuity issue across the interface. The accomplishment of the transmission conditions (in our case the field continuity) is ensured by the finite element contribution that does not compromise the required continuity.

Fig. 11 represents the functions F i (x) and G i (y) related to the separated representation contribution in domain Ω 1 . The ones concerning the domain Ω 2 can be obtained from these ones: functions of the x-coordinate are derived by symmetry considerations, and the ones concerning the y-coordinate are exactly the same that the corresponding ones in domain Ω 1 .

By combining the separated representation and the finite element contribution we can define the solution in the whole domain. The combined solution is depicted in Fig. 12. Fig. 13 depicts the error in the resulting solution for different number of terms in the finite sums decomposition of the separated representations (the finite element mesh remains unchanged). The solutions obtained by using a single term in the separated representations define the so-called Hartree approximations. Finally, Fig. 14 depicts the evolution of the L 2 -norm of the error as a function of the number of terms involved in the separated representation. As expected, the error stabilizes when the number of terms increases, and further reduction needs both finer separated representations (involving finer one-dimensional meshes) and finer finite element meshes.

Conclusions and perspectives

The proposed strategy can be used for enriching separated representations in order to take into account localized behavior of the solution. A numerical example proved the ability of defining reduced approximation bases able to capture localized behavior. This strategy can be also used for treating models involving strong or weak discontinuities. In this case the separated representation should be coupled with an appropriate finite element enrichment able to describe the involved discontinuities. There is an extensive literature concerning the use of finite elements for this purpose. There are many other potential applications. One could use this technique to define bridges between finite element models and other making use of proper generalized decompositions, or to make compatible subdomains making use of proper generalized decompositions involving different functional products or different one dimensional meshes on the common interface. It could be also applied for enforcing general boundary conditions by introducing a finite element description in a strip along the domain boundary, or to define proper generalized decompositions in geometrically complex domains. The analysis of these and other possibilities will be addressed in future works.

Figure 1 :

 1 Figure 1: Finite element mesh.

Figure 2 :

 2 Figure 2: Finite element solution.

Figure 3 :Figure 4 :

 34 Figure 3: Solution computed using a coarse separated representation.

Figure 6 :

 6 Figure 6: Computed solution using a finite element enriched separated representation.

Figure 7 :

 7 Figure 7: Finite element approximation contribution.

Figure 8 :

 8 Figure 8: Domain decomposition.

Figure 9 :

 9 Figure 9: Reference solution.

Figure 10 :

 10 Figure 10: Contributions to the whole solution: separated representation in Ω 1 (left); finite element contribution (center) and separated representation in Ω 2 (right).

Figure 11 :

 11 Figure 11: Functions F i (x) (left) and G i (y) involved by the separated representation in domain Ω 1 .

Figure 12 :

 12 Figure 12: Computed solution in the whole domain resulting of the combination of contributions to the solution depicted in Fig. 10.

Figure 13 :

 13 Figure 13: Error when the separated representations involve 1 (left), 3 (center) and 10 terms (right).

Figure 14 :

 14 Figure 14: Error versus number of terms involved in the separated representations.

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.


A W K Ma

F Chinesta

M R Mackley

The rheology and modeling of chemically treated carbon nanotubes suspensions

This paper reports recent experimental findings and rheological modeling on chemically treated single-walled carbon nanotubes ͑CNTs͒ suspended within an epoxy resin. When a CNT suspension was subject to a steady shear flow, it exhibited a shear-thinning characteristic, which was subsequently modeled by a Fokker-Planck ͑FP͒ based orientation model. The model assumes that the shear flow aligns CNT in the flow direction, but there are events such as Brownian motion and tube-tube interaction trying to randomize the orientation. In the FP orientation model, randomizing events were modeled with an appropriate rotary diffusion coefficient ͑D r ͒ and the shear-thinning behavior was explained in terms of progressive alignment of CNTs toward the shear direction. In terms of linear viscoelasticity ͑LVE͒, small-amplitude oscillatory measurements revealed mild elasticity for semidilute treated CNT suspensions. The exact origin for this elasticity is not clear and both tube-tube interaction and bending/stretching of CNTs have been proposed by other authors as possible origins. It is, however, clear from the current modeling that the experimental evolution of storage modulus ͑GЈ͒ cannot be described using a single-mode Maxwell model or simple Brownian rod modeling. In this paper, experimental LVE data of the treated CNT suspensions were fitted using the FP orientation model with an "effective diffusion coefficient" term and an empirical relation was subsequently identified for the effective diffusion term. Intuitively, chemical treatment has created a weakly interconnected network of CNT and it is believed that the mild elasticity originated from this weak network as well as other randomizing events ͑Brownian motion and tube-tube hydrodynamic interaction͒. Finally, step strain experiments confirmed the presence of a weak network at small strains, which at large strains was found to be destroyed. Incorporation of a strain softening factor allowed for the formulation of a 1 self-consistent FP based orientation model describing both the steady shear and LVE responses of treated CNT suspensions.

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 and nanodevices ͓Ajayan et al. 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.

III. STEADY SHEAR FLOW

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͔͒

= = -PI = + 2D = + = f = -PI = + =, ͑3͒
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͔͒,

dp គ dt = ⍀ = • p គ + kD = • p គ -k͑D =:͑p គ p គ ͒͒p គ , ͑4͒
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,

Ͷ S͑0,1͒ ͑x គ,p គ ,t͒dp គ = 1, ∀ x គ, ∀ t. ͑5͒
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,

a =͑x គ,t͒ = Ͷ S͑0,1͒ p គ p គ ͑x គ,p គ ,t͒dp គ ͑6͒
and the fourth order moment also known as the fourth order orientation tensor,

a = = ͑x គ,t͒ = Ͷ S͑0,1͒ p គ p គ p គ p គ ͑x គ,p គ ,t͒dp គ . ͑7͒
The evolution of the orientation distribution function is governed by the FP equation ͓see, for example, Hinch and Leal ͑1972͔͒,

d dt + ץ ץp គ • ͩ dp គ dt ͪ = ץ ץp គ • ͩ D r ץ ץp គ ͪ ,

͑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 dt = ץ ץt + v គ • grad͑͒, ͑9͒ where grad͑͒ = ץ ץx គ . ͑10͒
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

ץ ץt ͩͶ S͑0,1͒ ͑p គ ,t͒dp គ ͪ = 0, ∀ t,

͑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͔͒,

= f = 2N p ͑a = = :D =͒ + ␤D r ͩa = -I = 3 ͪ, ͑12͒
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

da = dt = ⍀ = • a = -a = • ⍀ = + k͑D = • a = + a = • D = -2͑a = = :D =͒͒ -6D r ͩa = -I = 3 ͪ. ͑13͒
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,

v គ = u v w = ␥ ˙y 0 0 .

͑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

a = 12 ␥ ˙= ͑1 + 2N p a 1212 ͒.

͑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,

D r0 = 3k B T͑ln͑L/d͒ -0.8͒ s L 3 , ͑17͒
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͒,

S = ͱ 3 2 ͑Q = :Q = ͒, ͑18͒
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͔͒.

IV. SMALL-AMPLITUDE LINEAR VISCOELASTIC OSCILLATORY DATA

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͒ and5͑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:

a ijkl lin = - 1 35
͑␦ ij ␦ kl + ␦ ik ␦ jl + ␦ il ␦ jk ͒ + 1 7 ͑a ij ␦ kl + a ik ␦ jl + a il ␦ jk + a kl ␦ ij + a jl ␦ ik + a jk ␦ il ͒.
͑a = = :D =͒ 12 Ϸ ␥ 15 . ͑20͒
Combining Eqs. ͑13͒ and ͑20͒ gives

da = 12 dt = ␥ 5 -6D r a = 12 , ͑21͒
which can be rearranged as da dt

+ a = ␥ ˙, ͑22͒
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 = 12 ͒ are

aЈ = ␥ 0 2 1 + 2 2 , aЉ = ␥ 0 1 + 2 2 .
͑23͒

The shear stress according to Eq. ͑12͒ is

Ј + iЉ = ␥ 0 iͩ1 + 2N p 15 ͪ+␤D r ͑aЈ + iaЉ͒ ͑ 24͒
and the storage modulus ͑GЈ͒ and the loss modulus ͑GЉ͒ are therefore

GЈ = ␤D r 2 1 + 2 2 , ͑25͒ GЉ = ͩ1 + 2N p 15 ͪ+␤D r 1 + 2 2 . ͑26͒
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

D ˆr = D r1 + ͑␥͒D r2 , ͑27͒
where

͑␥͒ = Ά ␥ c -␥ ␥ c if ␥ Յ ␥ c 0 if ␥ Ͼ ␥ c . · ͑28͒ D r1
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

= f = 2N p ͑a = = :D =͒ + ␤͑D r1 + ͑␥͒D r2 ͒ͩa = - I = 3 ͪ. ͑29͒
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 and8 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

da 12 dt = -6D r a 12 . ͑31͒
If D r is assumed to be independent of a 12 , integration of Eq. ͑31͒ would give

a 12 = Ce -6D r t , ͑32͒
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:

= f = 2N p ͑a = = :D =͒ + ␤͑D r1 + ͑␥͒D r2 ͒ͩa = - I = 3 ͪ, Fokker-Planck equation: d dt + ץ ץp គ ͩ dp គ dt ͪ = ץ ץp គ ͩ ͑D r1 + ͑␥͒D r2 ͒ ץ ץp គ ͪ ,
where

͑␥ ¯͒ = Ά ␥ c -␥ ␥ c if ␥ Յ ␥ c 0 if ␥ Ͼ ␥ c , · and
D r2 = C ͱ .
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.

L͑•͒

= d͑•͒ dt + ٌ sph • ͑p គ ͑•͒͒ -D r ⌬ sph ͑•͒ ͑ A2͒
and

ٌ sph • ͑p គ ͒ = 1 sin ץ ץ ͑sin ͓p គ ͔ ͒ + 1 sin ץ ץ ͓͑p គ ͔ ͒, ͑A3͒
where ͓p គ ͔ and ͓p គ ͔ refer to the and components of vector p គ , respectively; and finally

⌬ sph = 1 sin ץ ץ ͩsi n ץ ץ ͪ+ 1 sin 2 ץ 2 ץ 2 . ͑A4͒
If the finite difference method is applied and a grid on ⍀ is defined, consisting of N ϫ M nodes, then each grid node will be characterized by two indices ͑i , j͒ ͑1 Յ i Յ N and 1 Յ j Յ M͒ and the associated coordinates ͑ i , j ͒ will be given as i = ͑i -1͒h and j = ͑j -1͒h , where h =2 / N and h = / ͑M -1͒.

The coordinate transformation is given by x = sin cos , y = sin sin , ͑A5͒ z = cos .

Now, the discrete form of the FP equation is considered at each internal grid point ͑2 Յ i Յ N -1 and 2 Յ j Յ M -1͒, assuming the following finite difference formulas:

ͯ ͒•͑ץ ץ ͯ i,j = ͑•͉͒ i+1,j -͑•͉͒ i-1,j 2h , ͯ ͒•͑ץ ץ ͯ i,j = ͑•͉͒ i,j+1 -͑•͉͒ i,j-1 2h , ͑A6͒ ͯ ץ 2 ͑•͒ ץ 2 ͯ i,j = ͑•͉͒ i+1,j -2͑•͉͒ i,j + ͑•͉͒ i-1,j h 2 , ͯ ץ 2 ͑•͒ ץ 2 ͯ i,j = ͑•͉͒ i,j+1 -2͑•͉͒ i,j + ͑•͉͒ i,j-1 h 2 .
Time discretization is performed using an appropriate difference formula. If the nth time step ͑n Ն 1͒ is denoted by the superscript n ͓i.e., ͑•͒ i,j n = ͑•͒͑ = i , = j , t = n⌬t͔͒, the time derivative can be approximated by

ͯ ͒•͑ץ ץt ͯ i,j n = ͑•͒ i,j n -͑•͒ i,j n-1 ⌬t . ͑A7͒
At the boundary nodes consisting of i =1 or i = N finite difference formulas can also be defined taking into account the periodicity, that is, for ͑i , j͒ with i =1,

ͯ ͒•͑ץ ץ ͯ 1,j = ͑•͉͒ 2,j -͑•͉͒ N,j 2h , ͑A8͒ ͯ ץ 2 ͑•͒ ץ 2 ͯ 1,j = ͑•͉͒ 2,j -2͑•͉͒ 1,j + ͑•͉͒ N,j h 2 ,
and for ͑i , j͒ with i = N,

ͯ ͒•͑ץ ץ ͯ N,j = ͑•͉͒ 1,j -͑•͉͒ N-1,j 2h , ͑A9͒ ͯ ץ 2 ͑•͒ ץ 2 ͯ N,j = ͑•͉͒ 1,j -2͑•͉͒ N,j + ͑•͉͒ N-1,j h 2 .
Now, the finite difference form of the FP equation can be written for each node ͑i , j͒ with 1 Յ i Յ N and 2 Յ j Յ M -1, at each time t = n⌬t ͑n Ն 1͒, L͉͑͒ i,j n = 0 ∀ ͑i, j͒ ͓1, ... ,N͔ ϫ ͓2, ... ,M -1͔, ͑A10͒

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

͚ i=1 N ͚ j=1 M A ij i,j n = 1, ͑A12͒
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. 1 .FIG. 2 .

 12 FIG. 1. Optical micrograph of ͑a͒ 0.33% chemically treated and ͑b͒ 0.33% untreated CNTs suspended in epoxy resin ͑optical depth= 130 m͒.

  FIG. 3.Apparent viscosity ͑ a ͒ as a function of shear rate ͑␥ ˙͒ for ͑a͒ 0.33% chemically treated and ͑b͒ 0.33% untreated CNT suspensions.

  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 .

 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.

FIG. 7 .

 7 FIG. 7.Sensitivity of ͑a͒ GЈ, ͑b͒ GЉ, and ͑c͒ ء to fitting parameter C. ͑N p =10, ␤ =1000.͒

FIG. 8 .

 8 FIG.8. Sensitivity of ͑a͒ GЈ, ͑b͒ GЉ, and ͑c͒ ء to fitting parameter ␤. ͑N p =10, C = 0.05 s -0.5 . ͒

  FIG. 9. ͑a͒ Stress relaxation after a finite step strain ͑␥ 0 = 100%͒ was applied to different concentrations of CNT suspensions. ͑b͒ Stress relaxation data for the 0.5% CNT suspension with varying magnitudes of step strain.

=

  ∀ i ͓1, ... ,N͔.

  ͑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.

ACKNOWLEDGMENTS

The authors would like to thank Nanocomposites Inc. for providing the carbon nanotubes. A. Ma would also like to acknowledge the Croucher Foundation Scholarship and the Overseas Research Students Awards Scheme ͑ORSAS͒ for providing financial support.

APPENDIX: NUMERICAL ISSUES

To solve the FP equation, there are two simple options to discretize the unit sphere surface: ͑i͒ the use of Cartesian coordinates or ͑ii͒ the use of a spherical coordinates system. In the latter case, the equation is defined in the domain ͑ , ͒ ⍀ ϵ͓0,2͒ ϫ ͓0,͔. This essentially simplifies the application of discretization techniques, such as finite differences and finite elements, and allows for a computation with higher accuracy. The main difficulty associated with the use of spherical coordinates is the singularity at angles = 0 and = , but this is not a problem if it is defined that = 0 is perpendicular to both the shear and vorticity directions.

By assuming that the diffusion coefficient does not depend on the conformation coordinates, the spherical representation of the FP equation reads as