# Effects of a Bent Structure on the Linear Viscoelastic Response of Carbon Nanotube Diluted Suspensions

Camilo Cruz, Lounès Illoul, Francisco Chinesta, Gilles Regnier

► To cite this version:

Camilo Cruz, Lounès Illoul, Francisco Chinesta, Gilles Regnier. Effects of a Bent Structure on the Linear Viscoelastic Response of Carbon Nanotube Diluted Suspensions. Rheologica Acta, 2010, 49, pp.1141-1155. 10.1007/s00397-010-0487-0. hal-01006779

HAL Id: hal-01006779

https://hal.science/hal-01006779

Submitted on 18 Mar 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Effects of a bent structure on the linear viscoelastic response of diluted carbon nanotube suspensions

C. Cruz · L. Illoul · F. Chinesta · G. Régnier

**Abstract** Commonly isolated carbon nanotubes in suspension have been modelled as a perfectly straight structure. Nevertheless, single-wall carbon nanotubes (SWNTs) contain naturally side-wall defects and, in consequence, natural bent configurations. Hence, a semi-flexile filament model with a natural bent configuration was proposed to represent physically the SWNT structure. This continuous model was discretized as a non-freely jointed multi-bead-rod system with a natural bent configuration. Using a Brownian dynamics algorithm the dynamical mechanical contribution to the linear viscoelastic response of naturally bent SWNTs in dilute suspension was simulated. The dynamics of such system shows the apparition of new relaxation processes at intermediate frequencies characterized mainly by the activation of a mild elasticity. Storage modulus evolution at those intermediate frequencies strongly depends on the flexibility of the system, given by the rigidity constant of the bending potential and the number of constitutive rods.

**Keywords** Brownian dynamics · Semi-flexible filament · Linear viscoelasticity · Carbon nanotubes

## Introduction

Today, without any doubt nanotechnology is a hot topic. Belonging to this nano-world, carbon nanotubes are being intensively studied because their potential applications are extremely varied and the expected properties really impressive. A carbon nanotube (CNT) is a cylindrical structure constituted essentially by rolled-up graphite sheets with diameters going from 0.4 nm to several nanometres and lengths up to some millimetres (Pan et al. 1998). This particular structure confers to CNTs an extraordinary compromise between mechanical, thermal and electrical properties.

Particularly, single-walled carbon nanotubes (SWNTs) are considered as a perfect constituent material for obtaining very strong and ultra light electromechanical components based on nano-filaments (Chico et al. 1996; McEuen 2000), carbon fibres (Vigolo et al. 2000), and transparent films (Wu et al. 2004). Additionally, SWNTs have a great biomedical potential as an effective gene and drug delivery system through cell membranes (Shim et al. 2002; Kam et al. 2004).

Pristine SWNTs exhibit very strong tube-tube van der Waals forces that produce their natural agglomeration in bundles. In consequence, solubilisation of SWNTs as individual entities is difficult and remains the main obstacle to develop their industrial applications (Ajayan 1999). To solve this problem, chemical modification of the tube side wall by covalent bonding has been the most extended attempt to enhance the dispersion of SWNTs (Bahr and Tour 2001).

<footnote>C. Cruz · L. Illoul · G. Régnier (✉)
Laboratoire de Procédés et Ingénierie en Mécanique et Matériaux (PIMM, UMR CNRS 8006),
Arts et Métiers ParisTech, 151 bd de l'Hôpital,
75013 Paris, France
e-mail: gilles.regnier@ensam.eu</footnote>

<footnote>C. Cruz
e-mail: camilo.cruz@ensam.eu</footnote>

<footnote>F. Chinesta
EADS Corporate International Chair, Ecole Centrale
de Nantes, 1 rue de la Noë, BP 92101,
44321 Nantes Cedex 3, France</footnote>
A large range of the applications of SWNTs requires a transformation processing in liquid phase. For such reason, surface-treated SWNTs are usually suspended into a Newtonian fluid or a viscoelastic matrix (melted polymer). Control and optimization of those transforming processes, based on the flow of such suspensions, need a deep comprehension of their rheology (Kharchenko et al. 2004).

Davis et al. (2004) reported the dispersion of SWNTs in superacids (102% sulphuric acid), where the sidewall SWNT protonation is the driving mechanism that counteracts the tube-tube van der Waals forces and promotes the dispersion of individual SWNTs. They found that at dilute regime (concentration <120 ppm vol) SWNTs act as simple Brownian rods and the suspension exhibits an expected shear thinning behaviour.

In 2006, Rahatekar et al. studied the rheology of multi-wall carbon nanotubes (MWNTs) within an epoxy resin. They verified a two decades enhancement of the low shear viscosity when incorporating carbon nanotubes into the epoxy. When suspension was progressively sheared the suspension evidenced a rheo-thinning behaviour, a phenomenon mainly explained by the reduction of aggregates and not by the nanotube alignment. On the contrary, Fan and Advani (2005, 2007) revealed by transmission electron microscopy that MWNTs are aligned in the flow direction when their suspensions were sheared. They found also that the extent of such alignment is function of the shear rate. A more recent work showed that the shear-thinning behaviour of CNT suspensions can be explained considering both the aggregation and the alignment effects (Ma et al. 2008a).

In 2005, Xu et al. studied the morphology and rheology of carbon nanofiber suspensions. Carbon nanofibers (diameter of about 100–200 nm) are bigger than CNTs, but their structural similarities let to compare their dynamical behaviour. Xu et al. used several standard elastic and rigid dumbbell models to predict the rheological behaviour of aqueous carbon nanofiber suspensions. They conjectured that the elasticity in such suspensions comes mainly from the recoverable bending of the carbon nanofibers.

Recently, Ma et al. (2008b) presented a study about the extensional rheology of surface-treated SWNTs diluted within an epoxy matrix. They validated a filament thinning experiment to determine the extensional viscosity of CNT suspensions. It was observed that the filament relaxation process for the surface-treated SWNT samples (low concentration, ≤0.3% wt.) follows the typical Newtonian filament relaxation behaviour where filament diameter decreases linearly with time.

On the contrary, the extensional viscosity of the treated SWNT suspensions did not show the classical relationship between the zero-shear viscosity and the extensional viscosity for Newtonian liquids, giving evidence of an enhanced extensional viscosity, probably due to the high degree of alignment of treated SWNTs in the flow direction. This behaviour was modelled by a simple Fokker-Planck equation-based orientation model.

Given the high aspect ratio and tensile rigidity of carbon nanotubes the classical mesoscopic–macroscopic approaches for fibre suspensions (Hinch and Leal 1972; Petrie 1999; Sepehr et al. 2004; Rajabian et al. 2005) can be employed to describe the rheological behaviour of CNT suspensions. However, diffusion coefficients usually depend on empirical relationships (Folgar and Tucker 1984; Larson 1999) and lack of physical meaning.

Ma et al. (2009) reported a generalized Fokker-Planck-based model to predict the rheology of surface-treated SWNTs suspended within an epoxy matrix in the dilute to semi-dilute concentration regime. The classical shear-thinning behaviour when the suspension is subjected to a steady shear flow was explained in terms of the competition between two phenomena: SWNT orientation in the flow direction and SWNT random misalignment due to Brownian motion (thermal interaction with the solvent). The shear-thinning data were properly fitted by a Fokker-Planck orientation model tuning the rotary diffusion coefficient. On the other hand, when the surface-treated SWNT suspension was submitted to a small amplitude oscillatory deformation a non-negligible elastic response was identified in the tested frequencies. Figure 1 shows the viscoelastic contribution of functionalized SWNTs to the dynamic response of a 0.2 wt.% SWNT–epoxy resin suspension. A constant-slope evolution for the storage and loss modulus is observed in the tested frequency interval. This kind of rheological behaviour was found characteristic of the SWNT suspensions from the dilute to the semi-dilute concentration regime. Ma et al. (2009) fitted those linear viscoelastic data with the Fokker-Planck-based orientation model, but using an empirical relationship for the rotary diffusion coefficient that depends on the applied frequency. Even though there had been several hypotheses to explain this behaviour as, for example, existence of a weak CNT network (Hough et al. 2004) and electrostatic interaction between surface-treated SWNTs, the physical origins of this mild elasticity in the linear viscoelasticity spectra are not clear yet.

SWNTs can contain naturally topological defects (Lijima et al. 1992). Several studies have showed the
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Viscoelastic contribution of treated SWNTs to the dynamic response of a 0.2 wt.% SWNT-epoxy resin suspension. SWNT contribution to the storage and loss modulus are calculated as $G'_{CNT} = G'_{suspension} - G'_{epoxy}$ and $G''_{CNT} = G''_{suspension} - \omega\eta_{epoxy}$, respectively. Epoxy resin is supposed to be Newtonian ($G'_{epoxy} \sim 0$ and $\eta_{epoxy} = 10$ Pa s). Courtesy of Dr. Anson Ma (Rice University)</figcaption></figure>

existence of those defects by using Fresnel projection microscopy (Binh et al. 2000), scanning tunnelling microscopy (STM; Ouyang et al. 2001) and voltage-contrast scanning electron microscopy (VC-SEM; Vijayaraghavan et al. 2010). Most common topological defects in SWNTs are due to the existence of a non-hexagonal carbon ring (e.g. pentagon-heptagon defect) in middle of the classical honeycomb-like lattice structure of a CNT. Molecular dynamics (MD) is largely employed to identify the structural consequences of this kind of defects. For example, it has been demonstrated that a bent structure is thermodynamically more stable than a straight one once a localized topological defect appears in the SWNT growth (Xue et al. 2009). Additionally, it has been shown that a single pentagon-heptagon pair can bend the SWNT tubular structure forming a localized junction (Ajayan et al. 1998) with angles varying from 0° to 34° (Han et al. 1998; Wako et al. 2008) depending on the distance between the pentagon and heptagon pair defect (Lambin and Meunier 1999).

So, in order to contribute elucidating the physical phenomena behind the mild elasticity exhibited by the dynamics of individual SWNTs this paper investigates the influence of the naturally bent SWNT structure

on the linear viscoelastic response of their dilute suspensions. To do so, we proposed an equivalent mechanical model that represents an individual SWNT as a semi-flexible filament (or a non-freely multi-bead-rod model in a discretized version) with a natural bent configuration and we evaluated by Brownian dynamics the linear viscoelastic response of this physical model into a diluted suspension.

## Modelling

Persistent length ($L_p$) considerations are necessary when modelling high aspect ratio systems that exhibit a very important resistance to tangential extension or compression (Morse 1998; Shankar et al. 2002). In absence of external forces (e.g. a flow field), persistent length is a characteristic length that relates the bending rigidity of the system to the thermal energy around it. In practice, this parameter determines a length scale over which the high aspect ratio system presents significant curvature due to thermal interactions with the surroundings. In consequence, a system of total length $L$ where $L/L_p \ll 1$ can be considered as a simple rigid rod and a system where $L/L_p \gg 1$ can be considered as a flexible filament. Traditionally, individual SWNTs have been modelled as Brownian rigid rods (Davis et al. 2004; Duggal and Pasquali 2006; Mendes et al. 2008; Ma et al. 2009).

Nevertheless, there are experimental results that lead to think that treated SWNTs can be considered as other than rigid-rod systems. Using a near-infrared fluorescence technique, thermal SWNT-shape fluctuations has been observed into a diluted SWNT aqueous suspension (Tsybolski et al. 2005). Additionally, Fakhri et al. (2009) demonstrated that the observed bending dynamics corresponds accurately with a semi-flexible filament model ($L \sim L_p$), where $L_p$ scales with the $d^3$ (d, diameter of SWNT).

If we are considering now the action of an external flow field and the natural existence of a bent structure in SWNTs due to the presence of side-wall defects, the persistent length concept requires a more precise definition. In any case, those new considerations can only reduce the length scale of observable bending relative to the thermal classical persistent length and, hence, in this paper a semi-flexible filament is supposed to be a good equivalent micromechanical model describing the dynamics of an individual SWNT submitted to a homogeneous flow field.

In this context, we insist on considering that a SWNT contains naturally side-wall defects producing a
deviation from the perfectly straight structure; i.e. the
pristine SWNT has a slightly curved structure. This
assumption is supported on two main facts:

* Even if there is a large discussion about the influence of sample preparation on the observed structure of CNTs, there are numerous images taken by SEM, STM, atomic force microscopy (AFM) and transmission electron microscopy that reveal clearly a natural bent structure in some CNT samples (Shaffer et al. 1998; Clauss 1999; Ouyang et al. 2001; Kinloch et al. 2002; Loos et al. 2005; Lamprecht et al. 2009; Vijayaraghavan et al. 2010). For example, in Fig. 2 an image taken by AFM of surface-treated SWNTs protruding from a photo-cured acrylic surface is presented; the bent structure of the SWNTs is easily appreciable.

* The natural existence of topological defects that plays an important role in the mechanical and thermo-electrical properties of CNTs (Ruoff et al. 2003). A clear relationship between those CNT structural instabilities (pentagon/heptagon defect, for example) and the intrinsic curvature of CNT structures has been established (Clinard et al. 1994; Dimitrakopulos et al. 1997; Ajayan et al. 1998; Han et al. 1998; Clauss 1999; Lambin and Meunier 1999; Huang and Choi 2008; Wako et al. 2008; Xue et al. 2009).

So, in order to test the dynamic mechanical response
of this physical model, we decided to use an inertialess
Langevin dynamics approach, also known as Brownian

<figure><img src="image_3.png" /><figcaption>**Fig. 2** AFM image of treated SWNTs protruding from a photo-cured acrylic surface (image area corresponds to a 5-µm-side square). Courtesy of Dr. Anson Ma (Rice University), Dr. Loren Picco (Bristol University) and Pr. Mervyn Miles (Bristol University)</figcaption></figure>

dynamics. The semi-flexible filament was discretized in
a non-freely jointed multi-bead-rod model as observed
in Fig. 3. This multi-bead-rod model is composed of
*n* beads with positions **r**ᵢ, connected by *n* − 1 rods of
length *a*, where *a* approaches theoretically the persis-
tence length under flow conditions of the system. The
SWNT physical length is given by the total length of the
multi-bead-rod model (*n* − 1) *a*. Beads are considered
the centres of hydrodynamic resistance, as usual for
all the mechanical models in polymer kinetic theory
(Bird et al. 1987), and also mark the existence of a
bent topological defect. On the other hand, rods are
supposed infinitely rigid. We consider the existence of
an internal bending potential that mimics the SWNT
bending flexibility. Bending potential is given by the
action of a hypothetical torsion spring between each
pair of rods.

Concerning particularly the simulation of a dynamic
mechanical test on a diluted suspension of SWNTs,
several hypotheses are established:

* Dynamic mechanical solicitation occurs at constant temperature.

* Individual SWNTs are embedded in a Newtonian matrix of viscosity ηs.

* Suspension is supposed to be diluted, in consequence, no physical interaction between SWNTs is considered.

* Applied flow field is considered homogeneous, i.e. the rate-of-strain tensor is supposed constant within the suspension (Bird et al. 1987).

* No external forces (electrostatic, magnetic or gravitational) play a role in the dynamics of the suspension. Possible electrostatic interaction between surface-treatment groups is neglected because it

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Non-freely jointed multi-bead-rod model composed of *n* beads and *n* − 1 rods of length *a*. Bending potential between rods in the multi-bead-rod model is mimicked with a hypothetical torsion spring</figcaption></figure>
has been shown that those repulsion forces do not produce an increment of the storage modulus at high frequencies (Ma et al. 2008c).

## Kinematic formulation

Keeping in mind all previous hypotheses, the kinematic description of the multi-bead-rod system is given by the next n stochastic equations (n force balances for each bead):

$$ 
\begin{gathered}
\mathbf{0} = \mathbf{F}_i^d + \mathbf{F}_i^b + \mathbf{F}_i^\phi - \sum_j \mathbf{n}_{ij} \lambda_j \\
\text{for } i = 1 \dots n \text{ and } j = 1 \dots n-1 \qquad (1)
\end{gathered}
$$

where $\mathbf{F}_i^d$ is the hydrodynamic drag force acting on bead i, $\mathbf{F}_i^b$ is the Brownian force acting on bead i, $\mathbf{F}_i^\phi$ is the internal bending potential force acting on bead i and $\lambda_j$ is a constraint force associated to the rod connecting bead j to j + 1. Operator $\mathbf{n}_{ij}$ is defined as follows:

$$ \mathbf{n}_{ij} = \mathbf{u}_j (\delta_{i,j+1} - \delta_{i,j}) \qquad (2) $$

in which $\mathbf{u}_j = (\mathbf{r}_{j+1} - \mathbf{r}_j)/a$ is a unit vector collinear to the rod connecting beads j and j + 1. In order to complete correctly the kinematic description of the multi-bead-rod chain in a Cartesian coordinates system is necessary to impose a set of n - 1 constraints over the bead velocities:

$$ \dot{\mathbf{r}}_i \cdot \mathbf{u}_i - \dot{\mathbf{r}}_{i+1} \cdot \mathbf{u}_i = 0 \quad \text{for } i = 1 \dots n-1 \qquad (3) $$

where $\dot{\mathbf{r}}_i$ is the instantaneous velocity of bead i. In summary, using the n equations in (1) plus the n - 1 equations in (3) a linear system of equations describing the motion of the multi-bead-rod system can be obtained, where the unknown variables are n instantaneous bead velocities $\dot{\mathbf{r}}_i$ plus n - 1 constraint forces $\lambda_j$. In what follows a formal mathematical description for each term in Eq. 1 is developed and some important conditions are emphasized given the time-integration algorithm implemented in this work:

* Hydrodynamic drag force $\mathbf{F}_i^d$

This force describes the resistance experienced by the bead as it moves through the fluid. One of the simplest ways to express this force is given by the Stoke's law. This law considers that the hydrodynamic drag force is proportional to the difference between the fluid velocity and the averaged bead velocity. A more general expression takes into account a hydrodynamic interaction term that comes from the interactions inside

the same system as it is moving. Mathematically, this force writes:

$$ \mathbf{F}_i^d = \zeta ( [\mathbf{\kappa} \cdot \mathbf{r}_i] - \dot{\mathbf{r}}_i - \dot{\mathbf{r}}_i^* ) \qquad (4) $$

where $\zeta$ is a second-order friction tensor, $\mathbf{\kappa}$ is the gradient of the homogeneous velocity field acting in all the fluid and $\dot{\mathbf{r}}_i^*$ is a correction velocity that accounts for the variation in the local homogeneous flow field due to the motion of the other beads in the same system. In this work, an isotropic second-order friction tensor is supposed, reducing this parameter to a scalar friction coefficient $\zeta$ (i.e. $\zeta = \zeta \mathbf{I}$), and hydrodynamic interaction is neglected ($\dot{\mathbf{r}}_i^* = 0$) due to the small strain imposed and the mild bending rigidity of SWNTs. Previous suppositions are also necessary when implementing a simple one-step integration scheme (Öttinger 1996), as is the case of the current work.

* Brownian force $\mathbf{F}_i^b$

The Brownian forces account for the thermal effects of the solvent molecules over the SWNTs. Taking into account that the time scale of this thermal activity is quite small compared with the time scale of SWNT motion, Brownian forces are mathematically computed as a stochastic process. In unconstrained systems, the set of Brownian forces must to follow a normal distribution satisfying the next conditions:

$$ \langle \mathbf{F}_i^b (t) \rangle = 0 \qquad (5) $$

$$ \langle \mathbf{F}_i^b (t) \otimes \mathbf{F}_i^b (t + \Delta t) \rangle \approx \frac{2\zeta k_B T}{\Delta t} \mathbf{I} \qquad (6) $$

where $k_B$ is the Boltzmann constant, T is the absolute temperature and $\Delta t$ is the time step of the simulation. In Brownian dynamics, the Brownian force generation is a critical aspect, especially when the system contains constraints. A proof of this subtlety is the intense discussion about the simulations of coarse-grained models in the singular limit between rigid links and very stiff springs (Fixman 1978, 1986; Hinch 1994; Öttinger 1994, 1996; Grassia et al. 1995; Grassia and Hinch 1996). In the present modelling rods are considered intrinsically inextensible and they do not intend to mimic the effect of a very stiff spring. In order to satisfy the fluctuation dissipation theorem for systems containing physical constraints a geometrical projection is done over the unconstrained Brownian forces obtained using Eqs. 5 and 6. Brownian forces must to satisfy the next condition:

$$ 0 = \sum_i \mathbf{F}_i^b \cdot \mathbf{n}_{ij} \quad \text{for } j = 1 \dots n-1 \qquad (7) $$
A complete description of the algorithm used to satisfy Eq. 7 is given elsewhere (Montesi et al. 2005). Projection of Brownian forces is absolutely necessary when an anisotropic friction tensor is considered or second or higher order integration schemes are employed.

* Internal bending potential force $\mathbf{F}_i^\phi$

Based on the viscoelasticity theory of concentrated solutions of semi-flexible polymers (Morse 1998), Pasquali and Morse (2002) defined a discrete bending potential for a multi-bead-rod model perfectly straight in equilibrium:

$$ \phi = -\frac{K_b}{a} \sum_{i=2}^{n-1} \mathbf{u}_i \cdot \mathbf{u}_{i-1} \qquad (8) $$

where $\phi$ is the bending potential and $K_b$ is a bending rigidity constant. Using Eq. 8, a natural definition of the bending force for each bead can be obtained:

$$ \mathbf{F}_k^\phi = -\frac{\partial \phi}{\partial \mathbf{r}_k} = \frac{K_b}{a} \sum_{i=2}^{n-1} \frac{\partial (\mathbf{u}_i \cdot \mathbf{u}_{i-1})}{\partial \mathbf{r}_k} \qquad (9) $$

Appendix A shows that the total bending force on bead $k$ only depends on the local chain configuration and that this force can be interpreted as the sum of independent contributions coming from all the two-rod sub-sections containing bead $k$.

Considering now the main hypothesis mentioned in this paper, to take into account the naturally bent configuration of the multi-bead-rod model at equilibrium, the bending potential must to be redefined in the next way:

$$ \phi = -\frac{K_b}{a} \sum_{i=2}^{n-1} \mathbf{Z}_i \mathbf{u}_i \cdot \mathbf{u}_{i-1} = -\frac{K_b}{a} \sum_{i=2}^{n-1} \mathbf{u}_i' \cdot \mathbf{u}_{i-1} \qquad (10) $$

where $\mathbf{Z}_i$ is a linear operator that rotates vector $\mathbf{u}_i$ of $\theta_i^{eq}$ (interior angle between $\mathbf{u}_i$ and $\mathbf{u}_{i-1}$ at the equilibrium configuration) and $\mathbf{u}_i'$ is the rotated vector $\mathbf{Z}_i \mathbf{u}_i$. In an analogous way, the bending force expression in (9) takes the next form:

$$ \mathbf{F}_k^\phi = -\frac{\partial \phi}{\partial \mathbf{r}_k} = \frac{K_b}{a} \sum_{i=2}^{n-1} \frac{\partial (\mathbf{u}_i' \cdot \mathbf{u}_{i-1})}{\partial \mathbf{r}_k} \qquad (11) $$

Considering that the total bending force in Eq. 11 can also be interpreted as the sum of independent contributions coming from all the two-rod sub-sections containing bead $k$, as shown in Appendix A, the next algorithm is proposed to compute the total bending forces in coherence with the bent configuration:

1. Decomposition of the multi-bead-rod (n beads) system into n-2 sub-section systems of two consecutive rods as shown in Fig. 4 (note that sub-section p is composed by beads p, p+1 and p+2).

<figure><img src="image_5.png" /><figcaption>Fig. 4 Multi-bead-rod system of n beads (in black) decomposed in n-2 sub-sections p of two consecutive rods (in gray)</figcaption></figure>

2. Calculation of the non-zero contributions to the bending forces over all beads composing each sub-section p as follows:

* Bending force on bead p of sub-section p as:

$$ \mathbf{F}_{p,p}^{\phi} = -\frac{K_b}{a^2} (\mathbf{I} - \mathbf{u}_p \otimes \mathbf{u}_p) \cdot \mathbf{u}_{p+1}' \qquad (12) $$

* Bending force on bead p+2 of sub-section p as:

$$ \mathbf{F}_{p+2,p}^{\phi} = \mathbf{Z}_{p+1}^{-1} \cdot \left[ \frac{K_b}{a^2} \left( \mathbf{I} - \mathbf{u}_{p+1}' \otimes \mathbf{u}_{p+1}' \right) \cdot \mathbf{u}_p \right] \qquad (13) $$

Bending force on bead p+2 of sub-section p has been rotated back (operator $\mathbf{Z}_{p+1}^{-1}$) in order to correct the direction of the force, in coherence with the actual bent configuration.

* As the bending potential in each sub-section is locally independent of the adjacent sub-sections, it must to satisfy mechanical equilibrium and hence bending force on bead p+1 of sub-section p writes:

$$ \mathbf{F}_{p+1,p}^{\phi} = -\mathbf{F}_{p,p}^{\phi} - \mathbf{F}_{p+2,p}^{\phi} \qquad (14) $$

3. Finally, the total bending force on bead k is calculated as follows:

$$ \mathbf{F}_k^\phi = \sum_{p=1}^{n-2} \mathbf{F}_{k,p}^\phi \qquad (15) $$

Once the instantaneous velocity $\dot{\mathbf{r}}_i$ of the system is known (by solving the system constituted by the $2n-1$ equations defined in (1) and (3)), the system position
is updated using a one-step algorithm that conserves the rod lengths throughout the time step. For example, the numerical time integration for a bead–rod chain confined in the plan xy is obtained as follows:

1. Calculation of angular velocities $\omega_i = [00\omega_i]$ of rods by resolving adequately the over-defined system:

$$ \dot{\mathbf{r}}_{i+1} - \dot{\mathbf{r}}_i = \omega_i \otimes a\mathbf{u}_{i-1} \quad \text{for } i = 1 \dots n - 1 \quad (16) $$

2. Update of the position of bead 1 by:

$$ \mathbf{r}_{1}^{t+\Delta t} = \dot{\mathbf{r}}_{1}^{t} \cdot \Delta t \quad (17) $$

3. Update of the director vectors of rods by:

$$ \theta_{i}^{t+\Delta t} = \theta_{i}^{t} + \omega_{i} \cdot \Delta t \quad \text{for } i = 1 \dots n - 1 \quad (18) $$

4. Update of the rest of bead positions by resolving consecutively from $i = 2$ to $n$:

$$ \mathbf{r}_{i}^{t+\Delta t} = \mathbf{r}_{i-1}^{t+\Delta t} + a\mathbf{u}_{i-1}^{t+\Delta t} \quad (19) $$

Previous one-step algorithm is coherent with the explicit corrector–predictor scheme (Liu 1989; Somasi et al. 2002) because it ensures that the rotational velocity of each rod at the end of the time step is perpendicular to its orientation (Doyle et al. 1997).

## Stress tensor calculus

The total shear stress tensor $\tau$ in a suspension is supposed to be the sum of the contribution coming from the solvent $\tau_s$ and another one coming from the suspended entities $\tau_p$ (Bird et al. 1987):

$$ \tau = \tau_s + \tau_p \quad (20) $$

As defined in classical Newtonian fluid mechanics, $\tau_s$ in Eq. 20 can be also written as follows:

$$ \tau_s = \eta_s \dot{\mathbf{y}} \quad (21) $$

where $\eta_s$ is the solvent viscosity and $\dot{\mathbf{y}}$ is the homogeneous rate-of-strain tensor. To model the shear dynamic response, we are particularly interested on the shear stress contribution coming from the SWNTs (represented as non-freely multi-bead–rod chains) $\tau_p$. In order to quantify this contribution in the context of Brownian dynamics we used the Kramers–Kirkwood expression for the shear stress tensor:

$$ \tau_p = c \sum_v \langle \mathbf{R}_v \otimes \mathbf{F}_v^d \rangle \quad (22) $$

where c is the number of SWNTs per volume and $\mathbf{R}_v = \mathbf{r}_v - \mathbf{r}_c$ is the relative location of bead v to the centre of mass of the multi-bead–rod chain $\mathbf{r}_c$.

In order to correctly account for the Brownian stress contribution during a time step the shear stress tensor contribution coming from the SWNTs $\tau_p$ is calculated using the next algorithm:

1. Calculation of the hydrodynamic drag forces $\mathbf{F}_i^d$, the geometrical projected Brownian forces $\mathbf{F}_i^b$ and the internal bending potential forces $\mathbf{F}_i^\phi$ at time t as described in the kinematic description section.

2. Update of the system position at time $t + \Delta t$.

3. Using the same geometrically projected Brownian forces obtained at time t, recalculation of the system velocity and hydrodynamic drag forces $\mathbf{F}_i^{d,+\Delta t}$ at the updated configuration $t + \Delta t$.

4. Determination of the current shear stress tensor at time t by computing:

$$ \tau_p^t = c \sum_v \langle \mathbf{R}_v^{+\Delta t} \otimes \mathbf{F}_v^{d,+\Delta t} \rangle \quad (23) $$

## Computational details

Previous algorithm was programmed on C++, compiled in Linux 64-bits and run on a cluster 6 x Quad-Core AMD Opteron™ 2376 (2300 MHz) processors.

## Validation

In order to validate the Brownian dynamics algorithm we decided to simulate the dynamic response of the freely jointed three-bead model, also known as the trumbbell, trimer or (three-bead)–(two-rod) system. The dynamic response of bead–rod models for macromolecular solutions has been studied before by Hassager (1974). His results are in coherence with the Fokker-Plank equation for a bead–rod chain in kinetic theory (Kramers 1944; Bird et al. 1987), actually deviated from that one for bead–spring chains (Gottlieb and Bird 1976). Using asymptotic expansions, Hassager found a numerical solution for the trimer model, i.e. a system composed of three beads with friction coefficient $\zeta$ connected by two freely jointed rigid massless rods of length a. To simulate the trimer system using our Brownian dynamics algorithm we used a bidimensional modelling and we neglected the internal bending potential by establishing the bending rigidity constant $K_b$ equal to 0. To compare the Hassager’s results with our Brownian dynamics simulation, complex modulus is represented in a proper reduced frequency scale.

In Fig. 5 the Hassager’s numerical solution is compared with results of the previously described Brownian dynamics simulation for the case of a trimer system.
<figure><img src="image_6.png" /><figcaption>Fig. 5 The dynamic modulus of a freely jointed (three-bead)–(two-rod) system in a proper reduced frequency scale: comparison between the Hassager’s numerical solution and the current Brownian dynamics (BD) simulation. $τ_1$ is the main Hassager’s relaxation time and $λ_{rod}$ is the rotational diffusion time of an equivalent multi-bead rigid-rod system of equal length. BD parameters: $n = 2, c = 1, a = 1, ζ = 1, k_B T = 0.1$ and $K_b = 0$</figcaption></figure>

Hassager’s complex modulus is presented in a reduced frequency scale $ωτ_1$ where $τ_1$ is the main relaxation time given by the numerical solution (Hassager 1974):

$$τ_1 = \frac{\zeta a^2}{5.4376k_B T} \qquad (24)$$

On the other hand, our Brownian dynamics simulation for the trimer system is plotted in a reduced frequency scale $ωλ_{rod}$, where $λ_{rod}$ is the rotational diffusion time of a (three-bead)–(one-rod) system of equal length (Bird et al. 1987):

$$λ_{rod} = \frac{\zeta a^2}{3k_B T} \qquad (25)$$

The Brownian dynamics algorithm yielded a complex modulus in good agreement with the Hassager’s numerical solution for the trimer system. As bi-dimensional modelling validated the Hassager’s results, it was not necessary to implement the three-dimensional modelling, which would be more expensive in computational time. Figure 6 presents the performance of the Brownian dynamics code in function of the population considered. The performance of the Brownian dynamics code is measured in terms of the mean of relative errors for 17 storage modulus values homogeneously distributed in an interval of reduced frequencies going

<figure><img src="image_7.png" /><figcaption>Fig. 6 Mean relative error of the storage modulus (for a trimer system) and total computational time in function of the Brownian dynamics population. Mean relative errors are calculated as the average of 17 points distributed homogeneously in a reduced frequency range going from $10^{-2}$ to $10^2$. Relative error is calculated as follows: $\left|\left(\log G'_{simul}(\omega_R) - \log G'_{Hassager}(\omega_R)\right) / \log(G'_{Hassager}(\omega_R))\right|$</figcaption></figure>

from $10^{-2}$ to $10^2$. Loss modulus signal is practically not sensitive to the Brownian dynamics population; mean relative error for the loss modulus signal is 2.9%. A population of $10^4$ systems is selected for all the Brownian dynamics simulations because a good compromise between computational time and rheological prediction performance is guaranteed.

## Results and discussion

The Brownian dynamics algorithm described in the modelling section was used to predict the dynamic response of a non-freely jointed multi-bead-rod system with a natural bent configuration, a physical model intended to emulate the mechanical response under flow conditions of a SWNT that contains bendable junctions.

In the first stage, we decided to analyze the simplest case, i.e. the (three-bead)–(two-rod) system. The freely jointed version of this model (i.e. with $K_b = 0$) is presented in Fig. 7. In a general way, two different regions are clearly differentiated in the dynamic response of a freely jointed (three-bead)–(two-rod) system: a low frequency regime exhibiting a typical Maxwell terminal behaviour, where the thermal activity (or Brownian forces) drives the kinematic of the system; and a high frequency regime, where the domination of flow
<figure><img src="image_8.png" /><figcaption>Fig. 7 Predicted complex modulus responses of different (three-bead)–(two-rod) systems: freely jointed, non-freely jointed with natural straight configuration ($K_b = 0.1$ and $K_b = 10$) and non-freely-jointed with natural bent configuration ($K_b = 10$ and $\Delta\theta^{\max} = 30^{\circ}$). BD parameters: $n = 2$, $c = 1$, $a = 1$, $\zeta = 1$, $k_B T = 0.1$. $\lambda_{rod}$ is the rotational diffusion time of an equivalent multi-bead rigid-rod system of equal length</figcaption></figure>

forces over the diffusivity terms produces a plateau for the storage modulus and a limiting constant viscosity ($G''(\omega)/\omega$) – $\eta_s$ larger than zero. The existence of a non-zero limiting viscosity is physically more coherent than the classical Rouse-Zimm responses, where the model contribution to the dynamic viscosity disappears at high frequencies.

If we consider now the action of the bending potential, i.e. a non-freely jointed (three-bead)–(two-rod) system, a large range of dynamic responses is deployed in function of the equilibrium configuration and the bending rigidity constant. To explain this fact, in Fig. 7 the dynamic responses of a freely jointed system, two non-freely jointed systems with natural straight configuration and a non-freely jointed system with a naturally bent configuration are compared. Rod misalignment (i.e. the natural bent SWNT configuration) is generated following the next formula:

$$ \theta_i^0 = \theta_{i-1}^0 - \Delta\theta^{\max} + U(0, 1) \cdot 2\Delta\theta^{\max} \quad (26) $$

where $\theta_i^0$ is the equilibrium director vector of rod i, $\Delta\theta^{\max}$ is the maximal misalignment between two consecutive rods and $U(0, 1)$ is a continuous uniform random distribution function defined between 0 and 1. In a general way, when the bending rigidity constant $K_b$ is equal or greater than the thermal coefficient $k_B T$,

the loss modulus is enhanced over the entire frequency range and the storage modulus is increased in the low frequency regime. Otherwise, if the bending rigidity constant is much lower than the thermal coefficient ($K_b \ll k_B T$), a dynamic response similar to that one of a freely jointed system is obtained, showing just one dominant relaxation time associated with the thermal-to-flow-induced motion transition.

Additionally, in Fig. 7 the apparition of new relaxation processes when the bending rigidity constant is considerably greater than the thermal coefficient (i.e. $K_b \gg k_B T$) is observed. For instance, the activation of a mild elasticity at intermediate frequencies is clearly noticed in the dynamic response of the non-freely jointed (three-bead)–(two-rod) system with bending rigidity constant $K_b = 10$. When comparing the straight configuration with the bent one, even though the loss modulus responses are similar, an interesting difference in the storage modulus is appreciated. In fact, on the one hand the activation of the mild elasticity for the system with a natural bent configuration is observed at a lower characteristic frequency than the system with a straight configuration; on the other hand a more steeped mild elasticity is noticed for the bent system in comparison with the straight one. In front of this

<figure><img src="image_9.png" /><figcaption>Fig. 8 Predicted complex modulus responses of a non-freely jointed (three-bead)–(two-rod) system with a natural bent configuration (maximal misalignment angle $\Delta\theta^{\max} = 30^{\circ}$) for different values of the bending rigidity constant $K_b$. Freely jointed system ($K_b = 0$) is presented as reference. BD parameters: $n = 2$, $c = 1$, $a = 1$, $\zeta = 1$, $k_B T = 0.1$. $\lambda_{rod}$ is the rotational diffusion time of an equivalent multi-bead rigid-rod system of equal length</figcaption></figure>
qualitative feature the appropriateness of considering a
natural bent configuration in SWNTs is a priori justified
when analyzing its dynamic response under flow.

In Fig. 8 the Brownian dynamics predictions of
the complex modulus for a non-freely jointed
(three-bead)–(two-rod) system with a natural bent
configuration for different values of the bending
rigidity constant are presented. In Fig. 8 the typical
thermal-to-flow-induced relaxation process occurring
around the rotational diffusion time λ<sub>rod</sub> (as in the
equivalent straight systems) is easily identifiable.
Additionally, an interesting second relaxation process
(as showed before) activating a mild elasticity at
intermediate frequencies is again appreciated for the
natural bent systems. Brownian dynamics simulations
have showed that the characteristic frequency of this
relaxation process and the magnitude of the storage
modulus enhancement are function of the bending
rigidity constant. The activation frequency of this mild
elasticity is increased as the bending rigidity constant
is increased as well. This feature is physically coherent
taking into account that a higher energy input (i.e.
solicitation frequency) is required to activate the first
bending mode of a stiffer system. On the other hand, it
is also noticed that the storage modulus enhancement

<figure><img src="image_10.png" /><figcaption>**Fig. 10** Predicted loss modulus of a non-freely jointed multi-bead-rod system with a natural bent configuration (maximal misalignment angle $\Delta\theta^{\max} = 30^{\circ}$) for different number of rods $n$. BD parameters: $c = 1$, $a = 1$, $\zeta = 1$, $K_b = 1$, $k_B T = 0.1$. Loss modulus is scaled with the value of storage modulus at the unit reduced frequency. $\lambda_{rod}$ is the rotational diffusion time of an equivalent multi-bead rigid-rod system of equal length. *Verti-cal pointed line* represents the transition between the thermal-induced to the flow-induced regime</figcaption></figure>

and the frequency range associated to this new
relaxation process also rise as the bending rigidity
increases. This qualitative feature can be physically
explained considering that the energy stocking capacity
of a bending mode is directly proportional to the
bending rigidity constant.

The reduced storage modulus predictions of a non-
freely jointed multi-bead-rod system with a natural
bent configuration for different number of rods (n =
3, n = 5 and n = 10) are showed in Fig. 9. As in a
Kramers chain, in the present modelling the number
of rods reflects the flexibility of the bead-rod sys-
tem, but in addition it captures the density of bent
defects in the SWNT structure. In other words, this
parameter is directly related to the SWNT persistent
length under flow conditions that, at the same time,
we suggest strongly depends on the presence of bent
defects in the SWNT structure. For information, an ex-
ample of the equilibrium internal angles Δθ₀⁻ⱼ = θᵢ₀⁻
θⱼ₀ between consecutive rods used for the simulations
in Figs. 9 and 10 is presented in Table 1. Differences
in the dynamic response are only appreciable in the
flow-induced regime, i.e. at frequencies greater than
the characteristic relaxation frequency of an equivalent
multi-bead rigid-rod system. In other words, different

<figure><img src="image_11.png" /><figcaption>**Fig. 9** Predicted storage modulus of a non-freely jointed multi-bead-rod system with a natural bent configuration (maximal misalignment angle $\Delta\theta^{\max} = 30^{\circ}$) for different number of rods $n$. BD parameters are: $c = 1$, $a = 1$, $\zeta = 1$, $K_b = 1$, $k_B T = 0.1$. Storage modulus is scaled with the value of storage modulus at the unit reduced frequency. $\lambda_{rod}$ is the rotational diffusion time of an equivalent multi-bead rigid-rod system of equal length. *Verti-cal pointed line* represents the transition between the thermal-induced to the flow-induced regime</figcaption></figure>
<table><caption>**Table 1** Internal angles between consecutive rods ($\Delta\theta_{i-j}^0 = \theta_i^0 - \theta_j^0$) of the multi-bead-rod systems with natural bent configurations presented in Figs. 9 and 10</caption><thead><tr><th></th><th>Δθ<sup>0</sup><sub>2-1</sub></th><th>Δθ<sup>0</sup><sub>3-2</sub></th><th>Δθ<sup>0</sup><sub>4-3</sub></th><th>Δθ<sup>0</sup><sub>5-4</sub></th><th>Δθ<sup>0</sup><sub>6-5</sub></th><th>Δθ<sup>0</sup><sub>7-6</sub></th><th>Δθ<sup>0</sup><sub>8-7</sub></th><th>Δθ<sup>0</sup><sub>9-8</sub></th><th>Δθ<sup>0</sup><sub>10-9</sub></th></tr></thead><tbody><tr><td>3 rods</td><td>18.9</td><td>24.3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>5 rods</td><td>-22.4</td><td>24.8</td><td>7.9</td><td>-24.1</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>10 rods</td><td>-13.3</td><td>2.8</td><td>27.4</td><td>27.9</td><td>-20.5</td><td>28.2</td><td>27.4</td><td>-0.9</td><td>18.0</td></tr></tbody></table>

relaxation processes at intermediate frequencies are activated in function of the number of rods, leading invariably to different curves of the storage modulus. An enhancement of the elastic response when increasing the number of rods is observed and different slopes for the storage modulus in function of frequency are obtained (from 0.5 for $n = 3$ to 1.2 for $n = 10$ at reduced frequencies between $10^1$ and $10^2$). This mild elasticity can be explained by the coupled contribution of the tension forces required to maintain the inextensibility condition (as in the Kramers chain) and the forces coming from the torques activated by the bending of the system.

On the other hand, in Fig. 10 the reduced loss modulus responses for the same systems deployed on Fig. 9 are presented. In the low frequency regime the typical thermal motion response is observed for the three systems, i.e. a loss modulus evolving proportional with frequency. As in the storage modulus case, differences between the loss modulus responses for the three multi-bead-rod systems ($n = 3, n = 5$ and $n = 10$) are only appreciable in the high frequency range. Differences are quite less important than those found for the storage modulus because just light variations in the slopes were encountered (from 0.7 for $n = 3$ to 0.8 for $n = 10$ at reduced frequencies between $10^1$ and $10^2$).

In order to properly compare the experiments on Fig. 1 with our simulations we decided to represent the experimental data in a reduced scale. To do so, an estimation of the equivalent rotational diffusion time for the actual surface-treated SWNT is needed. The rotational diffusion time for a multi-bead rigid-rod defined by Bird et al. (1987) writes as follows:

$$ \lambda_{\text{rod}} = \frac{\zeta L^2 n(n+1)}{72(n-1)k_B T} \qquad (27) $$

where the isotropic bead friction coefficient $\zeta$ is a simple discretization of the bulk friction coefficient $\zeta_{\text{bulk}} = n\zeta$. The bulk friction coefficient for a rigid rod system of length $L$ and diameter $d$ can be calculated following the next equation (Mendes et al. 2008):

$$ \zeta_{\text{bulk}} = \frac{2\pi \eta_s L}{\log[L/d]} \qquad (28) $$

Rheological tests presented in Fig. 1 were effectuated at 25°C using a suspension of HiPco SWNT (Nanocomposites, Houston, TX, USA) within an epoxy resin with a Newtonian viscosity of 10 Pa s. Length of HiPco SWNTs largely depends on factors such as whether the tubes have been functionalized and how long they have been sonicated for. From the literature, we found that HiPco SWNTs have average lengths ranging from 400 to 700 nm and diameters between 0.8 and 1.3 nm (Parra-Vasquez et al. 2007; Mendes et al. 2008; Fakhri et al. 2009). Using mean values for the physical dimensions of HiPco SWNTs and employing Eqs. 27 and 28, we found equivalent rotational diffusion times for the SWNT going from 5.5 s for $n = 10$ to 9.0 s for $n = 3$.

In Fig. 11 the scaled experimental data are compared with the dynamical response of a naturally bent non-freely jointed multi-bead-rod system of five rods. It appears from Fig. 11 that the mild elasticity measured experimentally is a phenomenon occurring in the intermediate frequency range and that the simulation results produce coherent magnitudes. Observe, for

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Comparison in a proper reduced scale of the experimental SWNT dynamic response in diluted suspension (Fig. 1) with the BD predicted complex modulus of a naturally bent non-freely jointed multi-bead-rod system of five rods (Figs. 9 and 10)</figcaption></figure>
example, that the predicted ratio between G' and G'' at intermediate frequencies is similar to that one encountered in experimental dynamical tests of SWNT diluted suspensions. A more precise model fitting will require identifying the relaxation time spectra of an individual SWNT and the dynamic response of the diluted suspension in a higher frequency range. We emphasize that the principal objective of this paper is to demonstrate that a natural bent SWNT structure (due to presence of synthesis defects) can play an important role in the dynamical response of a SWNT diluted suspension and constitutes, in consequence, a physical parameter to take into account when modelling the mechanical behaviour of SWNTs.

The physical model proposed in this paper could be applied for simulating non-linear regimes, but probably we should consider intra- and inter-tube interactions, that can be introduced via the hydrodynamic interaction hypothesis (as done in polymer kinetic theory) coupled with an adapted Brownian dynamics strategy. The analysis of this extended model constitutes a work in progress. The results reported in this paper consider small amplitude oscillatory flows where the maximum deformation was chosen within the interval in which the results were independent on the maximum deformation applied, in order to ensure the validity of the linearity hypothesis.

## Conclusions

Isolated treated SWNTs have been modelled as semi-flexible filaments with a naturally bent configuration considering the probable existence of side-wall defects on carbon nanotubes. Continuous semi-flexible filament model was discretized as a non-freely jointed multi-bead-rod system with a non-straight configuration at equilibrium in order to implement a Brownian dynamics simulation. In this context, a mechanically coherent re-definition of the discretized version of the bending potential for a natural straight semi-flexible filament is required in order to tackle the cases of non-straight configurations.

New relaxation processes in the dynamic mechanical response of a naturally bent system are observed from the simulation results. Those processes are characterised by the activation of a mild elasticity at intermediate frequencies (revealed as an enhancement of the storage modulus). Moreover, non-negligible differences are found with respect to the dynamic response of a naturally straight system. It is observed that the enhancement of the storage modulus and the range of frequencies associated with those new relaxation

processes depend directly on the bending rigidity constant or, in other words, on the local structural stiffness of the SWNT junctions. Additionally, the dynamic response of such model is function of the persistent length under flow conditions (i.e. number of rods in the discretized version). An augmentation of the storage modulus slope in the intermediate frequency range is observed when increasing the number of rods.

This work will be useful to give richer explanations to the physical origins of the impressive stiffening effects coming from carbon nanotubes on the rheological behaviour of their diluted, semi-diluted and even concentrated suspensions. Future work should be pointed towards the identification of physical parameters and the consideration of more complex bending potentials and poly-dispersed length distributions.

**Acknowledgements** Authors would like to thank Dr. Anson Ma (Rice University) and Prof. Malcolm Mackley (University of Cambridge) for fruitful discussions and for providing helpful experimental information (Fig. 1). We would like to acknowledge Dr. Anson Ma (Rice University), Dr. Loren Picco (Bristol University) and Prof. Mervyn Miles (Bristol University) for allowing us to use the AFM image (Fig. 2). We would like to thank also the GeM institute (Institute de Recherche en Génie Civil et Mécanique) at Nantes (France) for providing the computational resources.

## Appendix A

Total bending force in a multi-bead-rod system writes:

$$ \mathbf{F}_{k}^{\phi}=-\frac{\partial \phi}{\partial \mathbf{r}_{k}}=\frac{K_{b}}{a} \sum_{i=2}^{n-1} \frac{\partial\left(\mathbf{u}_{i} \cdot \mathbf{u}_{i-1}\right)}{\partial \mathbf{r}_{k}} \quad (29) $$

Expanding Eq. 29 we have:

$$ \mathbf{F}_{k}^{\phi}=\frac{K_{b}}{a}\left(\frac{\partial(\mathbf{u}_{2} \cdot \mathbf{u}_{1})}{\partial \mathbf{r}_{k}}+\cdots+\frac{\partial(\mathbf{u}_{k-1} \cdot \mathbf{u}_{k-2})}{\partial \mathbf{r}_{k}}+\frac{\partial(\mathbf{u}_{k} \cdot \mathbf{u}_{k-1})}{\partial \mathbf{r}_{k}}+\frac{\partial(\mathbf{u}_{k+1} \cdot \mathbf{u}_{k})}{\partial \mathbf{r}_{k}}+\cdots+\frac{\partial(\mathbf{u}_{n-1} \cdot \mathbf{u}_{n-2})}{\partial \mathbf{r}_{k}}\right) \quad (30) $$

$$ \mathbf{F}_{k}^{\phi}=\frac{K_{b}}{a}\left(\frac{\partial \mathbf{u}_{2}}{\partial \mathbf{r}_{k}} \mathbf{u}_{1}+\frac{\partial \mathbf{u}_{1}}{\partial \mathbf{r}_{k}} \mathbf{u}_{2}+\cdots+\frac{\partial \mathbf{u}_{k-1}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k-2}+\frac{\partial \mathbf{u}_{k-2}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k-1}+\frac{\partial \mathbf{u}_{k}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k-1}+\frac{\partial \mathbf{u}_{k-1}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k}+\frac{\partial \mathbf{u}_{k+1}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k}+\frac{\partial \mathbf{u}_{k}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k+1}+\cdots+\frac{\partial \mathbf{u}_{n-1}}{\partial \mathbf{r}_{k}} \mathbf{u}_{n-2}+\frac{\partial \mathbf{u}_{n-2}}{\partial \mathbf{r}_{k}} \mathbf{u}_{n-1}\right) \quad (31) $$
Equation 31 can be evaluated using the identity:

$$ \frac{\partial}{\partial \mathbf{r}_k} \mathbf{u}_i = \frac{1}{a} (\delta_{k,i+1} - \delta_{k,i}) (\mathbf{I} - \mathbf{u}_i \otimes \mathbf{u}_i) \quad (32) $$

Observing identity (32) is clear that the derivative $\frac{\partial}{\partial \mathbf{r}_k} \mathbf{u}_i$ takes a non-zero value only when $i$ is $k-1$ or $k$. Hence, considering a bead $k$ in the middle of an infinite chain, Eq. 31 becomes:

$$ \mathbf{F}_{k}^{\phi} = \frac{K_{b}}{a} \left( \frac{\partial \mathbf{u}_{k-1}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k-2} \right) + \frac{K_{b}}{a} \left( \frac{\partial \mathbf{u}_{k}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k-1} + \frac{\partial \mathbf{u}_{k-1}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k} \right) \\ + \frac{K_{b}}{a} \left( \frac{\partial \mathbf{u}_{k}}{\partial \mathbf{r}_{k}} \mathbf{u}_{k+1} \right) \quad (33) $$

Total bending force on bead $k$, as given in the previous equation, can be interpreted as the sum of three independent contributions. In order to explain the origin of those contributions consider a multi-bead-rod system containing $n$ beads decomposed into $n-2$ independent sub-sections of two consecutive rods as showed in Fig. 4. Taking each of those sub-sections as an independent system, the expressions for the bending forces associated to each sub-section are given by Eq. 29. For example, sub-section $p$ (composed by beads $p, p+1$ and $p+2$) has the associated bending forces $\mathbf{F}_{p,p}^{\phi}$, $\mathbf{F}_{p+1,p}^{\phi}$ and $\mathbf{F}_{p+2,p}^{\phi}$, where first sub-index refers to the bead and second sub-index refers to the sub-section. Mathematically, those bending forces write:

$$ \mathbf{F}_{p,p}^{\phi} = \frac{K_b}{a} \frac{\partial (\mathbf{u}_{p+1} \cdot \mathbf{u}_p)}{\partial \mathbf{r}_p} = \frac{K_b}{a} \frac{\partial \mathbf{u}_p}{\partial \mathbf{r}_p} \mathbf{u}_{p+1} \quad (34) $$

$$ \mathbf{F}_{p+1,p}^{\phi} = \frac{K_b}{a} \frac{\partial (\mathbf{u}_{p+1} \cdot \mathbf{u}_p)}{\partial \mathbf{r}_{p+1}} = \frac{K_b}{a} \frac{\partial \mathbf{u}_{p+1}}{\partial \mathbf{r}_{p+1}} \mathbf{u}_p \\ + \frac{K_b}{a} \frac{\partial \mathbf{u}_p}{\partial \mathbf{r}_{p+1}} \mathbf{u}_{p+1} \quad (35) $$

$$ \mathbf{F}_{p+2,p}^{\phi} = \frac{K_b}{a} \frac{\partial (\mathbf{u}_{p+1} \cdot \mathbf{u}_p)}{\partial \mathbf{r}_{p+2}} = \frac{K_b}{a} \frac{\partial \mathbf{u}_{p+1}}{\partial \mathbf{r}_{p+2}} \mathbf{u}_p \quad (36) $$

Applying identity (32) to the three previous equations we have:

$$ \mathbf{F}_{p,p}^{\phi} = -\frac{K_b}{a^2} (\mathbf{I} - \mathbf{u}_p \otimes \mathbf{u}_p) \mathbf{u}_{p+1} \quad (37) $$

$$ \mathbf{F}_{p+1,p}^{\phi} = -\frac{K_b}{a^2} (\mathbf{I} - \mathbf{u}_{p+1} \otimes \mathbf{u}_{p+1}) \mathbf{u}_p \\ + \frac{K_b}{a^2} (\mathbf{I} - \mathbf{u}_p \otimes \mathbf{u}_p) \mathbf{u}_{p+1} \quad (38) $$

$$ \mathbf{F}_{p+2,p}^{\phi} = \frac{K_b}{a^2} (\mathbf{I} - \mathbf{u}_{p+1} \otimes \mathbf{u}_{p+1}) \mathbf{u}_p \quad (39) $$

From the previous equations it results clear that bending forces are in mechanical equilibrium. For this reason:

$$ \mathbf{F}_{p+1,p}^{\phi} = -\mathbf{F}_{p,p}^{\phi} - \mathbf{F}_{p+2,p}^{\phi} \quad (40) $$

Now, using the previous results about the decomposition of a multi-bead-rod system is easy to explain the origin of different contributions in Eq. 33 as follows:

$$ \mathbf{F}_{k}^{\phi} = \mathbf{F}_{k,k-2}^{\phi} + \mathbf{F}_{k,k-1}^{\phi} + \mathbf{F}_{k,k}^{\phi} \quad (41) $$

In other words, total bending force on bead $k$ can be interpreted as the sum of all the independent bending forces coming from sub-sections containing bead $k$.

## References

* Ajayan PM (1999) Nanotubes from carbon. Chem Rev 99(7): 1787-1799

* Ajayan PM, Ravikumar V, Charlier J (1998) Surface reconstruction and dimensional changes in single-walled carbon nanotubes. Phys Rev Lett 81(7):1437-1440

* Bahr JL, Tour JM (2001) Highly functionalized carbon nanotubes using in situ generated diazonium compounds. Chem Mater 13(11):3823-3824

* Binh VT, Vincent P, Feschet F, Bonard J (2000) Local analysis of the morphological properties of single-wall carbon nanotubes by Fresnel projection microscopy. J Appl Phys 88(6):3385-3391

* Bird RB, Curtiss CF, Armstrong RC, Hassager O (1987) Dynamics of polymer liquids. Vol. 2: Kinetic Theory. Wiley InterScience, New York

* Chico L, Crespi VH, Benedict LX, Louie SG, Cohen ML (1996) Pure carbon nanoscale devices: nanotube heterojunctions. Phys Rev Lett 76(6):971-974

* Clauss W (1999) Scanning tunneling microscopy of carbon nanotubes. Appl Phys A-Mater 69(3):275-281

* Clinard C, Rouzaud JN, Delpeux S, Beguin F, Conard J (1994) Electron microscopy, growth and defects of carbon nanotubes. J Phys Chem Solids 55(7):651-657

* Davis VA, Ericson LM, Parra-Vasquez NG, Fan H, Wang Y, Prieto V, Longoria JA, Ramesh S, Saini RK, Kittrel C, Billups WE, Adams WW, Hauge RH, Smalley RE, Pasquali M (2004) Phase behaviour and rheology of SWNTs in superacids. Macromolecules 37:154-160

* Dimitrakopulos GP, Dravid DP, Karakostas TH, Pond RC (1997) The defect character of carbon nanotubes and nanoparticles. Acta Crystallogr A 53(3):341-351
Doyle PS, Shaqfeh ESG, Gast AP (1997) Dynamic simulation of freely draining flexible polymers in steady linear flows. J Fluid Mech 334:251-291

Duggal R, Pasquali M (2006) Dynamics of individual single-walled carbon nanotubes in water by real-time visualization. Phys Rev Lett 96:246104

Fakhri N, Tsyboulski DA, Cognet L, Weisman RB, Pasquali M (2009) Diameter-dependent bending dynamics of single-walled carbon nanotubes in liquids. Proc Natl Acad Sci USA 106(34):14219-14223

Fan Z, Advani SG (2005) Characterization of orientation state of carbon nanotubes in shear flow. Polymer 46(14):5232-5240

Fan Z, Advani SG (2007) Rheology of multiwall carbon nanotube suspensions. J Rheol 51(4):585-604

Fixman M (1978) Simulation of polymer dynamics. I. General theory. J Chem Phys 69(4):1527-1537

Fixman M (1986) Implicit algorithm for Brownian dynamics of polymers. Macromolecules 19:1195-1204

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

Gottlieb M, Bird RB (1976) A molecular dynamics calculation to confirm the incorrectness of the random-walk distribution for describing the Kramers freely jointed bead-rod chain. J Chem Phys 65(6):2467-2468

Grassia P, Hinch EJ (1996) Computer simulations of polymer chain relaxation via Brownian motion. J Fluid Mech 308:255-288.

Grassia PS, Hinch EJ, Nitsche LC (1995) Computer simulations of Brownian motion of complex systems. J Fluid Mech 282:373-288

Han J, Anantram MP, Jaffe RL, Kong J, Dai H (1998) Observation and modelling of single-wall carbon nanotube junctions. Phys Rev B 57(23):14983-14989

Hassager O (1974) Kinetic theory and rheology of bead-rod models for macromolecular solutions. II. Linear unsteady flow properties. J Chem Phys 60(10):4001-4008

Hinch EJ (1994) Brownian motion with stiff bonds and rigid constraints. J Fluid Mech 271:219-234

Hinch EJ, Leal LG (1972) The effect of Brownian motion on the rheological properties of a suspension of nonspherical particles. J Fluid Mech 52:683-712

Hough LA, Islam MF, Janmey PA, Yodh AG (2004) Viscoelasticity of single wall carbon nanotube suspensions. Phys Rev Lett 93:168102

Huang J, Choi WB (2008) Controlled growth and electrical characterization of bent single-walled carbon nanotubes. Nanotechnology 19:505601

Kam NWS, Jessop TC, Wender PA, Dai H (2004) Nanotube molecular transporter: internalization of carbon nanotube-protein conjugates into mammalian cells. J Am Chem Soc 126(22):6850-6851

Kharchenko SB, Douglas JF, Obrzut J, Grulke EA, Migler KB (2004) Flow-induced properties of nanotube-filled polymer materials. Nat Mater 3:564-568

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

Kramers HA (1944) Het gedrag van macromoleculen in een stroomende vloeistof. Physica 11(1):1-19

Lambin PH, Meunier V (1999) Structural properties of junctions between two carbon nanotubes. Appl Phys A 68:263-266

Lamprecht C, Danzberger J, Lukanov P, Tilmaciu CM, Galibert AM, Soula B, Flahaut E, Gruber HJ, Hinterdorfer P, Ebner A, Kienberg F (2009) AMF imaging of function-

alized double-walled carbon nanotubes. Ultramicroscopy 109(8):899-906

Larson RG (1999) The structure and rheology of complex fluids. Oxford University Press, New York

Lijima S, Ichihashi T, Ando Y (1992) Pentagons, heptagons and negative curvature in graphite microtubule growth. Nature 356(6372):776-778

Liu TW (1989) Flexible polymer chain dynamics and rheological properties in steady flows. J Chem Phys 90:5826-5842

Loos J, Alexeev A, Grossiord N, Koning CE, Regev O (2005) Visualization of single-wall carbon nanotubes (SWNT) networks in conductive polystyrene nanocomposites by charge contrast imaging. Ultramicroscopy 104:160-167

Ma WKA, Chinesta F, Ammar A, Mackley MR (2008a) Rheological modeling of carbon nanotube aggregate suspensions. J Rheol 52(6):1311-1330

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

Ma A, Cruz C, Giner A, Mackley M, Régnier G, Chinesta F (2008c) Modelling elastic behaviour in functionalized carbon nanotubes suspensions. Inter J Mater Form, Suppl 1: 631-634

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

McEuen PL (2000) Single-wall carbon nanotubes. Phys World 46:1804-1811

Mendes MJ, Schmidt HK, Pasquali M (2008) Brownian dynamics simulations of single-wall carbon nanotube separation by type using electrophoresis. J Phys Chem B 112:7467-7477

Montesi A, Morse DC, Pasquali M (2005) Brownian dynamics algorithm for bead-rod semiflexible chain with anisotropic friction. J Chem Phys 122:084903

Morse DC (1998) Viscoelasticity of concentrated isotropic solutions of semiflexible polymers. 1. Model and stress tensor. Macromolecules 31(20):7030-7043

Öttinger HC (1994) Brownian dynamics of rigid polymer chains with hydrodynamic interactions. Phys Rev E 50(4): 2696-2701

Öttinger HC (1996) Stochastic processes in polymeric fluids. Springer, Berlin.

Ouyang M, Huang J, Cheung CL, Lieber CM (2001) Atomically resolved single-walled carbon nanotube intramolecular junctions. Science 291(5501):97-100

Pan ZW, Xie SS, Chang BH, Wang CY, Lu L, Liu W, Zhou WY, Li WZ, Qian LX (1998) Very long carbon nanotubes. Nature 394(6694):631-632

Parra-Vasquez ANG, Stepanek I, Davis VA, Moore VC, Haroz EH, Shaver J, Hauge RH, Smalley RE, Pasquali M (2007) Simple length determination of single-walled carbon nanotubes by viscosity measurements in dilute suspensions. Macromolecules 40(11):4043-4047

Pasquali M, Morse DC (2002) An efficient algorithm for metric correction forces in simulations of linear polymers with constrained bond lengths. J Chem Phys 116(5):1834-1838

Petrie CJS (1999) The rheology of fibre suspensions. J Non-Newton Fluid 87(2-3):369-402

Rahatekar SS, Koziol KKK, Butler SA, Elliott JA, Shaffer MSP, Mackley MR, Windle AH (2006) Optical microstructure and viscosity enhancement for an epoxy resin matrix containing multiwall carbon nanotubes. J Rheol 50(5):599-610

Rajabian M, Dubois C, Grmela M (2005) Suspensions of semiflexible fibers in polymeric fluids: rheology and thermodynamics. Rheol Acta 44(5):521-535
Ruoff RS, Qian D, Liu WK (2003) Mechanical properties of carbon nanotubes: theoretical predictions and experimental measurements. CR Phys 4:993–1008

Sepehr M, Carreau PJ, Grmela M, Ausias G, Lafleur PG (2004) Comparison of rheological properties of fiber suspensions with model predictions. J Polym Eng 24(6):579–610

Shaffer MSP, Fan X, Windle AH (1998) Dispersion and packing of carbon nanotubes. Carbon 36(11):1603–1612

Shankar V, Pasquali M, Morse DC (2002) Theory of linear viscoelasticity of semiflexible rods in dilute solution. J Rheol 46(5):1111–1154

Shim M, Kam NWS, Chen RJ, Li Y, Dai H (2002) Functionalization of carbon nanotubes for biocompatibility and biomolecular recognition. Nano Lett 2(4):285–288

Somasi M, Khomami B, Woo NJ, Hur JS, Shaqfeh ESG (2002) Brownian dynamics simulations of bead–rod and bead–spring chains: numerical algorithms and coarse-graining issues. J Non-Newton Fluid Mech 108(1–3):227–255

Tsybolski DA, Bachilo SM, Weisman RB (2005) Versatilisation on individual single-walled carbon nanotubes with near-infrared fluorescence microscopy. Nano Lett 5(5): 975–979

Vigolo B, Penicaud A, Coulon C, Sauder C, Pailler R, Journet C, Bernier P, Poulin P (2000) Macroscopic fibers and ribbons of oriented carbon nanotubes. Science 290(5495):1331–1334

Vijayaraghavan A, Marquardt CW, Dehm S, Hennrich F, Krupke R (2010) Imaging defects and junctions in single-walled carbon nanotubes by voltage-contrast scanning electron microscopy. Carbon 48:494–500

Wako K, Oda T, Tachibana M, Kojima K (2008) Bending deformation of single-walled carbon nanotubes caused by 5–7 pair couple defect. Jpn J Appl Phys 1 47(8):6601–6605

Wu Z, Chen Z, Du X, Logan JM, Sippel J, Nikolou M, Kamaras K, Reynolds JR, Tanner DB, Hebard AF, Rinzler AG (2004) Transparent, conductive carbon nanotube films. Science 305(5688):1273–1276

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

Xue B, Shao X, Cai W (2009) Comparison of the properties of bent and straight single-walled carbon nanotube intramolecular junctions. J Chem Theory Comput 5(6): 1554–1559