Chapter 13: Multiscale modeling of polymer–carbon nanotube composites – Polymer-Carbon Nanotube Composites


Multiscale modeling of polymer–carbon nanotube composites

G.M. Odegard,     Michigan Technological University USA


Carbon nanotube (CNT) composite materials have the potential to provide significant increases in specific stiffness and specific strength relative to materials used for many engineering structural applications. To facilitate the design and development of CNT composites, structure–property relationships must be established that predict the bulk thermo-mechanical response of these materials as a function of the molecular and micro-structure. The objective of this chapter is to describe a general framework for multiscale modeling of CNT composites. First, the fundamental aspects of efficient and accurate modeling techniques are discussed. This is followed by a review of current state-of-the-art modeling approaches. Finally, a specific example is presented that describes the application of the approach to a specific CNT composite system.

Key words

continuum mechanics


molecular dynamics

representative volume element

structure–property relationships

13.1 Introduction

While many experimental-based methods for characterizing the mechanical behavior of composite materials have proven to be robust and accurate, they incur high costs due to expensive testing equipment, specimen fabrication, and labor. Fortunately, computational modeling can be used to facilitate the development, characterization, and testing of composite materials by providing reliable and efficient structure–property relationships that are determined via inexpensive computational simulations.

In the past two decades, computational molecular modeling approaches (Leach, 2001) have emerged as important tools that can be used to predict atomic structure, vibrational frequencies, binding energies, heats of reaction, electrical properties, and mechanical properties of organic and inorganic materials. These methods are ideal for studying the behavior of atoms on the scale of nanometers or below. While the very largest of these simulations can include up to a billion (1×109) atoms, most engineering structures contain atoms numbering on the order of at least Avagadro’s number (1×1023 atoms). Clearly, many orders of magnitude exist between the number of atoms that computational chemistry methods can model and which engineering structures contain.

Since the 17th century, mathematicians, scientists, and engineers have continually worked to establish models for the understanding of behavior of bulk quantities of homogeneous and composite materials (Truesdell and Toupin, 1960; Herakovich, 1998; Truesdell and Noll, 2004). These models are important for understanding the motion, deformation, failure, buckling, and vibration of engineering structures and components. Because models that are established for these purposes cannot realistically contain all of the details used in computational chemistry simulations, approximations to the behavior of material must be made for which the number of atoms is minimally on the order of 1 × 1023 atoms. The assumption in continuum mechanics is that matter is modeled in a three-dimensional Euclidean point space for which points in a body in motion do not represent individual atoms. Instead, a body in motion contains an infinite number of points which form a continuum.

For the development of carbon nanotube (CNT) nanocomposite materials for engineering applications, it is often necessary to establish structure–property relationships, for which the molecular structure must be related to engineering-scale behavior. This clearly requires tools to relate the predicted molecular structure and behavior from computational chemistry techniques to the engineering structural level. These tools must span the length scale difference between Angstroms and meters. Many studies have focused on the modeling of CNT composites at more than one length-scale, including molecular and continuum-based scales (Frankland et al., 2003; Odegard et al., 2003; Odegard et al., 2004; Odegard et al., 2005b; Shi et al., 2005; Jiang et al., 2006; Grujicic et al., 2007; Karakasidis and Charitidis, 2007; Mokashi et al., 2007; Liu et al., 2008; Lu et al., 2008; Tserpes et al., 2008; Zeng et al., 2008; Shen, 2009; Takeda et al., 2009). These studies employ a wide range of philosophical approaches. A generalized framework needs to be established that provides a common structure to these different approaches.

This chapter presents generalized methods for multi-scale modeling of CNT composites. Basic descriptions of the computational chemistry and continuum mechanics modeling tools are described first. Following this, methods for linking the different length-scale models are presented. This chapter is not intended to present an exhaustive review of the modeling techniques of CNT composites that numerous researchers have reported in the past few years. Instead, this chapter focuses on establishing a mathematical framework that presents a general approach to modeling CNT/polymer composite materials.

13.2 Computational modeling tools

Molecular dynamics (MD) (Leach, 2001) is a computational technique in which a time evolution of a set of interacting atoms is followed by integrating their equations of motion. The forces between atoms are due to the interactions with the other atoms. A trajectory is calculated in a 6-N dimensional phase space (three position and three momentum components for each of the N atoms). Typical MD simulations of CNT composites are performed on molecular systems containing up to tens of thousands of atoms and for simulation times up to nanoseconds.

The physical quantities of the system are represented by averages over configurations distributed according to the chosen statistical ensemble. A trajectory obtained with MD provides such a set of configurations. Therefore the computation of a physical quantity is obtained as an arithmetic average of the instantaneous values. Statistical mechanics is the link between the nanometer behavior and thermodynamics. Thus the atomic system is expected to behave differently for different pressures and temperatures.

The interactions of the particular atom types are described by the total potential energy of the system, U, as a function of the positions of the individual atoms at a particular instant in time


where Xi represents the coordinates of atom i in a system of N atoms. The potential equation is invariant to the coordinate transformations, and is expressed in terms of the relative positions of the atoms with respect to each other, rather than from absolute coordinates.

Continuum mechanics (Truesdell and Toupin, 1960; Truesdell and Noll, 2004) describes the motion and interaction of a set of 0-dimensional particles that form a mathematical continuum. For a given set of initial and boundary conditions of a particular volume of a continuous medium, the motion and interaction of the particles are governed by a series of fundamental laws. The balance of mass governs the change in mass that occurs for a control volume. The balance of linear momentum and angular momentum for a continuum ultimately leads to equations of motion. The balance of energy forces the continuous system to follow the first law of thermodynamics and includes changes in energy due to damage or other physical phenomena. The Clausius-Duhem inequality ensures that the continuous system obeys the second law of thermodynamics. Finally, the constitutive equations govern the behavior between thermodynamic forces (e.g. temperature, applied deformation, applied electric and magnetic fields) and thermodynamics fluxes (e.g. stress, entropy, electric polarization).

Perhaps the most important aspect of continuum mechanics is its absolute flexibility for any type of material system. By appropriately modifying the balance laws and constitutive equations, complex material systems such as biological tissue, cellular materials, piezoelectric materials, and composites can be modeled on any length-scale level. Although this flexibility allows the behavior of any material to be modeled on the bulk-level, it also renders continuum models to be dependent on empirical data obtained through testing for accurate material characterization. This is in contrast to the MD approach in which the atomic behavior is established from basic physical principles of atomic theory.

Often, structures must be analyzed that are assumed to be composed of a continuous material that behaves according to a predefined continuum theory. For many complex structure geometries and loading conditions, closed-form solutions to the governing laws of continuum mechanics are not available. In this case, the finite element analysis (FEA) method (Bathe, 1996) is particularly useful in predicting the behavior for a wide range of structural geometries and load conditions. Using FEA, the stress and strain fields in an engineering structure are determined by discretizing the continuum into elements of primitive shapes (e.g. bricks, tetrahedrons). The nodes are at the corners, and sometimes, on the mid-sides of the element boundaries. As long as the geometry of the elements (mesh) is not too coarse, the overall predicted properties of FEA models can be accurately estimated by solving for the stress and strain fields of all of the elements in the model using standard numerical approaches.

Micromechanical methods are used to predict the bulk properties of continuous yet heterogeneous solids. While the assumption of continuity is maintained, it is assumed that different phases of the material exist that interact with each other by transferring load. A large number of micromechanical approaches have been developed (Herakovich, 1998) with a wide range of assumptions. While many approaches are efficient with less accuracy, others are mathematically complex with tremendous accuracy. It is important to note that FEA can be used as a micromechanical tool.

The analysis tools described in this section will be placed in the context of multiscale modeling throughout the rest of this chapter.

13.3 Equivalent-continuum modeling concepts

Many material properties are used to describe the behavior of continuum-based systems, such as Young’s modulus, and ultimate strength. While these types of parameters are indispensable to design engineers, they are defined for a mathematical continuum. Therefore, these quantities cannot be used to describe the properties of a discrete structure. For example, a lattice structure (truss, frame, or molecular structure) as a whole cannot have a corresponding value of Young’s modulus associated with it directly. However, the lattice structure can be modeled as an equivalent-continuum structure in which the overall behavior is similar to that of the lattice structure for the same loading conditions. The equivalent continuum must have equivalent material properties such as Young’s modulus associated with it in order to describe the mechanical behavior. This raises the question: How are the properties of an equivalent-continuum of a discrete structure determined? This is the topic of this section.

13.3.1 Representative volume element

Consider Fig. 13.1, which shows three different molecular models with different length-scale levels. The structure of these three molecular models is compared with the structure of an element of a mathematical continuum. When comparing the molecular model with a side length of 10 Å to the continuum model, it is clear that the points in the continuum do not have a one-to-one correspondence with the atoms in the molecular model. Even if the centroid of each atom in the molecular model were mapped to the corresponding points in the continuum model, there are points in the continuum model that do not correspond to atoms in the molecular model, thus eliminating a one-to-one correspondence of points and atoms. As the length scale of the molecular model increases from 10 Å to 50 Å, it is clear that there are more points in the continuum model that correspond to centroids of atoms in the molecular model, although there is still not a one-to-one correspondence between the two. As the length scale continues to grow to 200 Å, it becomes difficult to discern individual atoms in the image, as the atom sizes become very small compared to the molecular model size. In fact, the cubes of the molecular models with side lengths of 10 Å and 200 Å contain approximately 80 and 640,000 atoms, respectively. Therefore, comparison of the 200 Å molecular model with the continuum model clearly shows the nearly continuous nature of the 200 Å molecular model, as the significance of the behavior of individual atoms significantly decreases with increased molecular model size. In fact, it is generally assumed that molecular models of about the size of the 200 Å model behave like the continuous, mathematical continuum model with clearly defined elastic constants and smooth fields. However, large molecular models are computationally difficult to work with, so methods are necessary to model smaller molecular models as continuum-based entities.

13.1 Different length scales of molecular models (Odegard, 2009). (Reproduced with kind permission of Springer Science and Business Media.)

A representative volume element (RVE) is such a model. The morphological structural features of the RVE do not necessarily behave like points on a homogeneous continuum when a motion or deformation is applied. Instead, the RVE is defined as the volume element that efficiently and effectively describes the structure of the material in a statistical sense. That is, if a solid is built up of identical RVEs that are placed side-by-side with a perfect packing, then the resulting macrostructure will effectively describe the structure of the material. In some cases, RVEs can perfectly describe the minimum building block necessary to describe the continuum, such as the unit cell of a single crystal metal (Fig. 13.2) or the fiber composite material with hexagonal packing (Fig. 13.3). These materials are referred to as periodic materials. In other cases, the RVE can only approximately describe the minimum-sized building block of a material, such the RVE of bulk metallic glass (Fig. 13.4). Such amorphous materials have no long-range order, so that an RVE of an amorphous material can never exactly model the structure of the material at any spatial point. Rather, the RVE of an amorphous material only models the structure in a statistically accurate manner.

13.2 Crystal unit cell.

13.3 Fiber composite microstructure (Odegard, 2009). (Reproduced with kind permission of Springer Science and Business Media.)

13.4 Bulk metallic glass molecular model (Odegard, 2009). (Reproduced with kind permission of Springer Science and Business Media.)

In Figs 13.2, 13.3, and 13.4, the periodicity vectors are indicated by L(1), …, L(4). The periodicity vector components describe the length between opposite sides of a RVE. For the cubic RVEs shown in Figs 13.2 and 13.4, there are three vectors, L(1), L(2) and L(3), each corresponding to a pair of sides of the RVE. For the more complex RVE shown in Fig. 13.3, there are four vectors L(1), …, L(4).

The question remains: How is the RVE of a heterogeneous non-continuous structure used to establish the continuum-based balance laws and constitutive equations? As we have seen, the balance laws and constitutive equations of continuum mechanics assume the presence of a mathematical continuum, which clearly does not exist in the RVEs shown in Figs 13.2, 13.3, and 13.4. The answer to this question will be addressed in the next section.

13.3.2 Equivalent continuum

At some level, no material is truly described by a mathematical continuum. All materials contain a heterogeneous or discrete structure on some length-scale level. Some materials contain a definite structure at the micrometer length scale, such as wood and fiber-reinforced composites. Some materials exhibit a nearly homogeneous texture except on the atomic level, such as single-crystal metals. Therefore, it can be argued that any material that is described by a constitutive equation is an equivalent continuum. However, the concept of the equivalent continuum is applied to the modeling of materials at length scales in which the characteristic size of the inhomogeneity exists. For example, the modeling of the mechanical behavior of a fiber-reinforced composite material at length scales that are orders of magnitude larger than the fiber diameters does not require an effective continuum, as the material on that level closely approximates a continuum (recall the large molecular model of Fig. 13.1). That is, the fibers are so small with respect to the modeled continuum that they behave similar to a material point in the field equations described above. If the behavior is modeled at length scales that are about the same size of the fiber, then the material will no longer behave in the smooth manner associated with the field equations of a classical continuum. In this case, the equivalent-continuum properties are determined by studying the effective mechanical response of the RVE. Consider the molecular RVE shown in Fig. 13.5. Suppose that an equivalent volume and shape of a mathematically continuous solid is used to represent the RVE, denoted as region (with boundary ∂). The equivalent-continuum solid should mimic the behavior of the molecular model as closely as possible under all mechanical loadings. It is important to note that an equivalent-continuum model will ideally represent the behavior of an arbitrary volume of the actual material, not just the shape defined by the RVE. The definition of is only necessary to relate the mechanical response of the RVE to that of the equivalent material points in the equivalent continuum.

13.5 Equivalent-continuum model of a molecular RVE.

The mechanical response of the RVE can only be studied by applying boundary conditions to the RVE. This is usually performed computationally. For example, displacements can be applied to the surface of a RVE to calculate the resulting loads, and hence mechanical properties, of a solid material. Although electric and magnetic fields can also be applied as boundary conditions, the following discussion will focus on mechanical properties. There are four types of boundary conditions for RVEs in static equilibrium: (1) displacement-controlled boundary conditions (also called kinemetic or Dirichlet boundary conditions); (2) traction-controlled boundary conditions (also called static or Neumann boundary conditions); (3) periodic boundary conditions; or (4) a mixture of these three. The boundary conditions are applied at locations in the RVE (such as the centroids of atoms) that have the same locations as material points in the equivalent continuum .

For displacement-controlled boundary conditions, the components of the prescribed displacement vector ū (boldface indicates a vector quantity) are specified everywhere on ∂


where u is the displacement vector. For traction-controlled boundary conditions, the components of the prescribed traction vector are specified everywhere on ∂ as


where Sij are the stress tensor components and Nj are the surface unit normal vector components. For periodic boundary conditions, the prescribed displacements and tractions are given by


where uave is the average displacement vector associated with the bulk deformation of the solid material, and L is the periodicity vector of the RVE (shown in Figs 13.2, 13.3, and 13.4). For the case of mixed boundary conditions, the boundary of region can be divided into three sub-boundaries ∂d, ∂1t, and ∂p, such that


where Ø is the null set and the superscript ° denotes the relative interior. The corresponding boundary conditions are


Although the simplest approach to applying boundary conditions to a RVE (either computationally or analytically) is to use the displacement- or traction-controlled boundary conditions, it has been shown (Jiang et al., 2001; Jiang et al., 2002) that the application of periodic or mixed boundary conditions to the RVE results in a more realistic material response.

When determining effective-continuum properties of a material by computationally simulating the deformation of a RVE, the natural question arises: What size do the edges of a RVE need to be? The answer to this question depends mostly on two issues: The periodicity of the specific material and the boundary conditions that will be used to compute the material structure and properties. Specifically, the minimum size of an RVE for a periodic material is the minimum size of a possible repeatable structure, as shown in Fig. 13.3 for the smaller RVE of the fiber-composite material. However, if periodic boundary conditions are applied to the edges of the RVE, then the minimal size of the RVE is the minimum size necessary to construct the material structure without rigid rotations of the RVE, such as the larger RVE in Fig. 13.3 and the RVE in Fig. 13.2. If periodic boundary conditions are not used, then the computational results will depend on the RVE size. This dependence will depend on the morphology and properties of the material. Generally, the larger the RVE, the more accuracy is obtained with the results. For example, consider the elastic properties of aligned-fiber composites established by Jiang et al. (Jiang et al., 2002). Traction-controlled, displacement-controlled, and periodic boundary conditions were applied for bulk-level transverse shear (shear in the plane transverse to the fiber direction). Two cases were simulated, one which contained inclusions that were stiffer than the matrix by an order of magnitude and one in which the matrix was stiffer than the inclusions by an order of magnitude. Young’s modulus of the matrix was assumed to be 1 GPa. The corresponding transverse shear modulus was determined using all three boundary conditions for the two composite systems for RVE sizes of δ = δ0, 2δ0,3δ0 and 4δ0, where δ0 is the RVE size shown in the inset of Figs 13.6 and 13.7. Figures 13.6 and 13.7 show the calculated shear modulus for the stiff inclusions and matrix, respectively, for the different RVE lengths. The data has been plotted with smooth lines to emphasize the overall trend. Clearly, for the case of stiff inclusions shown in Fig. 13.6, the discrepancy between the predicted modulus from the periodic and displacement boundary conditions is larger than that between the periodic and traction boundary conditions. Both discrepancies decrease as the RVE size increases. For the case of the stiff matrix shown in Fig. 13.7, the predicted shear modulus from the displacement boundary conditions has better agreement with the periodic boundary conditions than does the modulus predicted with the traction boundary conditions. Again, as the RVE size increases, the different sets of boundary conditions predict a more similar shear modulus.

13.6 Transverse shear modulus of fiber composite with stiff inclusions (Odegard, 2009). BC, boundary condition. (Reproduced with kind permission of Springer Science and Business Media.)

13.7 Transverse shear modulus of fiber composite with stiff matrix (Odegard, 2009). (Reproduced with kind permission of Springer Science and Business Media.)

For amorphous materials, such as that shown in Fig. 13.4, the minimum size of the RVE is the size that statistically represents the structure of the material. This is a fairly ambiguous statement; however, there is no general agreement of a minimum necessary size of an RVE that produces more accurately predicted bulk-level properties. Of course, as the RVE size increases, the more likely the predicted properties will agree with those measured at the bulk-length scale level.

Once an RVE of a molecular structure is established, then it must be used to established equivalent-continuum constitutive properties. This process is described in the next section.

13.3.3 Equivalence of averaged scalar fields

So far, the necessity of constructing an RVE to mimic the bulk-scale mechanical response (equivalent continuum) of a material has been discussed. However, the determination of the properties of the equivalent continuum based upon the response of the RVE to applied boundary conditions has not been detailed. The establishment of an equivalent continuum model that behaves the same as the RVE under identical boundary conditions is important when attempting to understand molecular structure/bulk-level property relationships.

In general, an equivalent continuum must meet two requirements in order to accurately predict the behavior of a particular material:

1. Under identical applied far-field deformations (or loads), the RVE and the equivalent continuum must have identical (or nearly identical) values of one or more scalar fields that are averaged over the volume of the RVE and volume .

2. The material points of the equivalent continuum volume must have the same kinematic motion as material points (atoms in some cases) of the heterogeneous RVE at the same locations relative to some defined basis set and origin.

This subsection describes Requirement 1, the Equivalence of averaged scalar fields, in detail.

Requirement 1 expresses the need to have one or more scalar parameters; such as the scalar strain-energy density or the six scalar components of the stress tensor averaged over the RVE; to have equal values under identical loads applied to the RVE and equivalent-continuum models. For example, if periodic boundary conditions are applied to the RVE via Equation 13.4 such that the total potential energy, as calculated with an atomic potential, is equal to the strain energy of the effective continuum if a homogeneous deformation field, which matches the deformations averaged over the RVE, is applied to the effective continuum. In other words, if identical deformations are applied to a very large RVE and the equivalent continuum, the total energies should be the same. Moreover, the material parameters of the equivalent continuum can be adjusted such that the energies of the two models match. It is important to note that many researchers choose to match components of virial and continuum stress tenors (from the RVE and effective continuum models, respectively) under identical conditions. Although this approach makes intuitive sense, matching a single scalar parameter, such as strain energy, is a much more efficient approach than matching six independent components of a symmetric stress tensor. Further details on this requirement can be found elsewhere (Odegard et al., 2002). Once Requirement 1 is satisfied, then Requirement 2 must be considered.

13.3.4 Kinematic equivalence

Requirement 2 is often referred to as the Cauchy-Born Rule, which requires the kinematic motions of the RVE and the equivalent continuum to match for each atom (in the case of a molecular RVE) or each material point (for a heterogeneous micro-scale model) of the RVE. For example, in the case of the RVE shown in Fig. 13.2, the equivalent continuum should deform in the same identical manner on both the edges and in the interior as the center of the atoms in the RVE. This requirement often requires higher-order elasticity theories (Aifantis, 1984; Eringen, 1999, 2002) to be used to accurately match RVE and effective continuum deformations.

Although this requirement can be easily satisfied for the deformations of simple RVEs such as that shown in Fig 13.2, this rule is usually ignored for more complex RVEs, such as those shown in Figs 13.3 and 13.4, for one of the following reasons. Either (a) this requirement is over-restrictive and unnecessary given the required predictive accuracy of the effective continuum, or (b) this requirement is extremely difficult to impose on very complex RVEs. An example of point (a) is captured with the heterogeneous RVE of Fig. 13.3. In the composites community, it is often unnecessary to have predictive effective continuum models to predict the point-to-point kinematic mechanical behavior of the microstructure. The relaxation of this requirement has presented few difficulties in the successful design and implementation of most fiber-reinforced composite materials in the past several decades. An example of item (b) for the amorphous RVE is shown in Fig. 13.4. Given the complex atomic interactions that occur on this length scale for a set of atoms that have no local geometric order, the kinematic motion of the atoms is not expected to be uniform. Establishing a higher-order effective continuum model to match a highly non-uniform deformation field would be an exhausting and unnecessary task. Therefore, although Requirement 2 is rigorously followed for simple crystalline (or highly ordered) material systems, it is rarely followed for amorphous or structurally complex materials. The following section describes how these requirements are met for specific material systems.

13.4 Specific equivalent-continuum modeling methods

This section addresses specific approaches to establishing an equivalent-continuum representation of a molecular representation of a material.

13.4.1 Fluctuation methods

A broad class of methods has been developed to predict effective mechanical properties of materials based on the thermal fluctuation behavior of a RVE when no loads or deformations are applied. The resulting fluctuation methods are based on the work of Parrinello and Rahman (Parrinello and Rahman, 1982). As the RVE is simulated using MD, the geometric fluctuations of the RVE are correlated with the corresponding potential energy of the system. Because the average stress components in the RVE can be related to the derivative of free energy with respect to deformation, a series of stress–strain behaviors are observed computationally as the molecular system evolves over time. The thermal fluctuations of the RVE include axial deformation and shear deformation modes. Therefore, for a given sampling of simulation time, the six independent components of stress can be correlated with the six independent components of strain, and the corresponding elastic constants are subsequently established via the constitutive equation. A wide range of specific techniques have been utilized for specific materials and simulation ensembles (Ray and Rahman, 1984; Gusev et al., 1996; Zhou, 2001; Van Workum and de Pablo, 2003; Meyers et al., 2005; Yoshimoto et al., 2005; Miranda et al., 2006).

Although the fluctuation methods offer an efficient approach to establishing elastic constants of a wide range of materials, the obvious drawback is the inability of this approach to predict the mechanical response of materials to large deformations. To achieve this, the RVE must be subjected to applied loads or deformations of the appropriate magnitude. This process is described in the next two sub-sections for static and dynamic problems.

13.4.2 Static deformation methods

Perhaps the most common approach to determining the elastic properties of an effective continuum is with the static deformation method (Theodorou and Suter, 1986; Fan and Hsu, 1992; Hutnik et al., 1993; Fan et al., 1994; Raaska et al., 1994; Aleman and Munoz-Guerra, 1996; Kang et al., 1998; Jang and Jo, 1999a; Jang and Jo, 1999b; Yang and Jo, 2001; Odegard et al., 2003; Saether et al., 2003; Wang et al., 2003; Odegard et al., 2004; Odegard et al., 2005a; Odegard et al., 2005b; Li and Chou, 2006; Raffaini et al., 2006; Wang et al., 2006; Wu and Xu, 2006; Papakonstantopoulos et al., 2007; Valavala et al., 2007). Although many variations exist, this technique involves the application of homogeneous axial, volumetric and shear deformations/loads to the boundaries of an RVE. Consistent with requirement #1 of an effective continuum, the corresponding strain energy, stress, or strain is calculated. Constitutive equations are used to determine the elastic properties based on the response of the RVE to applied loads. This approach is applicable to amorphous material systems.

A reference configuration of the molecular structure of the RVE is determined by subjecting an equilibrated molecular structure to an energy minimization. The reference configuration of the RVE is subsequently exposed to the applied deformations. Initially, the boundaries and the atoms are displaced uniformly and the molecular structure is again subjected to an energy minimization to establish the deformed configuration. The potential energy or the stresses are calculated and compared to the corresponding values in the reference configuration, and elastic properties are determined based on the constitutive equations. Alternatively, loads can be applied to the reference configuration, and the resulting deformations can be determined for an equilibrated deformed state. This method has been extended into a hyperelastic framework to allow for large deformations of materials (Valavala et al., 2007).

The entire static molecular deformation process is usually performed using molecular mechanics or FEA. Therefore, the effects of the thermal motion of atoms on the mechanical behavior of the material are disregarded with this method. However, the efficiency and accuracy of this method have made it a popular choice for predicting bulk mechanical properties of finite-sized molecular systems. For modeling time-dependent effects, such as thermal motion, many researchers use the dynamic methods described in the next sub-section.

13.4.3 Dynamic deformation methods

Similar to the static deformation method, the dynamic deformation method involves deformation of the RVE to determine effective continuum properties (Raaska et al., 1994; Frankland et al., 2002; Brown et al., 2003; Frankland et al., 2003; Jin and Yuan, 2003; Griebel and Hamaekers, 2004; Sheng et al., 2004; Qi et al., 2005; Valavala et al., 2007). However, in the dynamic approach, the motion of the atoms in response to a deformation is determined with MD. Therefore, the time-dependent response of the atoms to applied deformations is computed with the integrated equations of motion, not energy minimizations. This allows the effect of the thermal motion of the atoms on the predicted properties to be determined.

Associated with the deformation of the RVE with this approach is a strain rate (or loading rate). The applied deformation is incrementally prescribed onto the RVE to simulate an applied bulk strain rate. However, because such simulated strain rates are orders of magnitude larger than those experienced in the laboratory and in most engineering applications, the resulting calculated effective continuum properties may not match experimentally-determined values and they may not be practical for the design of engineering materials. Furthermore, determination of effective continuum properties with the dynamic approach generally requires significantly more computational time than with the static approach. The selection of one of the methods described in this sub-section strongly depends on the material and behavior type that need to be modeled. Every particular problem has a different approach, as demonstrated in the following section.

13.5 Example: polymer–carbon nanotube composite

In 2005, Odegard et al. (Odegard et al., 2005b) reported on the constitutive modeling of CNT–polyethylene composites for the two cases in which the CNT reinforcement was functionalized and was not functionalized with the surrounding polyethylene molecules. The purpose of the investigation was to determine if the chemical bonds between the CNTs and surround polymer matrix transferred more load than the van der Waals forces alone in a non-functionalized composite. A SWNT (10,10) of radius 6.78 Å was modeled as being surrounded by polyethylene polymer chains. Both functionalized and non-functionalized systems were modeled. Figure 13.8 shows the modeled RVE for the two systems. In the composite with CNT functionalization, two polymer chains were attached to six carbons on the CNT by chemical linkages of 2 CH2 groups.

13.8 RVEs of non-functionalized and functionalized CNT composite materials (Odegard, 2009). e1, e2, e3 = basis set reactors. (Reproduced with kind permission of Springer Science and Business Media.)

The geometry of the equivalent continuum was assumed to be cylindrical for ease in use in subsequent micromechanical analyses. Thus, the equivalent continuum for these materials is henceforth referred to as the effective fiber. It was assumed that the polymer molecules that were near the polymer–CNT interface were included in the effective fiber. The effective fiber radius and length were 1.1 and 4.3 nm, respectively. The material composing the effective fiber was assumed to exhibit orthotropic symmetry. The elastic properties of the effective fiber were determined by equating the total strain energies of the effective fiber and molecular model RVE under identical boundary conditions. Further details on this portion of the analysis can be found elsewhere (Odegard et al., 2005b). Once the effective fiber properties were determined, a micromechanics analysis was performed to model the reinforcement of a polyethylene matrix reinforced with effective fibers (Fig. 13.8). It was assumed that the matrix polymer surrounding the effective fiber had mechanical properties equal to those of bulk amorphous polyethylene. All relative movement and interaction of the polymer chains with respect to each other were modeled at the molecular level. This motion and interactions were therefore indirectly considered in the subsequent determination of the properties for the effective fibers, and it is therefore assumed that perfect bonding existed between the CNT–polymer effective fibers and the surrounding polymer matrix in the micromechanical analysis.

The elastic stiffness components, volume fraction, length, and orientation of the effective fiber were used for the fiber properties in the micromechanical analysis. Although the CNT and effective-fiber lengths were equivalent, the CNT volume fraction was determined to be 52.9% of the effective-fiber volume fraction. The overall composite stiffness was calculated for effective fiber lengths up to 500 nm and effective fiber volume fractions up to 90%, which corresponds to the maximum volume fraction for hexagonally packed fibers. The calculations were performed assuming both perfectly aligned and three-dimensional randomly oriented effective fibers in an amorphous polyethylene matrix.

The calculated Young’s moduli of the CNT composites are plotted in Fig. 13.9 as a function of CNT length, for a 2% CNT volume fraction. The longitudinal Young’s modulus of the aligned composite, E1, and Young’s modulus of the random composite, E (E = E1 = E2 = E3), had a nonlinear dependence on the CNT length. An increase in Young’s modulus with respect to an increase in CNT length is expected to correspond to a relative increase in load transfer efficiency between the CNT and surrounding polymer. The data in Fig. 13.9 indicate that at a CNT length of about 250 nm, the efficiency of load transfer is nearly maximized, as evidenced by the nearly zero slope in the data curve. Further increases in CNT length beyond 250 nm resulted in relatively small increases in Young’s modulus for a given CNT volume fraction. As the CNT length became greater than approximately 100 nm, Young’s modulus for the composite without CNT functionalization became larger than that of the composite with CNT functionalization for the random composite and longitudinal deformation of the aligned composite. At 500 nm, the functionalization reduced the longitudinal Young’s modulus of the aligned composite and Young’s modulus of the random composite by 11 and 9%, respectively. In contrast, the transverse Young’s modulus of the aligned composite, E2 ≈ E3, had no dependence on CNT length. Also, there was no effect of functionalization on the transverse deformation of the aligned composite.

13.9 Young’s moduli of the composite systems for a 2% CNT volume fraction.

The longitudinal Young’s moduli of the random and aligned composites are plotted in Fig. 13.10 as a function of CNT volume fraction, for a constant CNT length of 500 nm. At this CNT length, the change in the longitudinal Young’s modulus with respect to CNT volume fraction of the aligned composite was nearly linear. Over the complete range of CNT volume fractions, the functionalization of the CNT reduced the longitudinal Young’s modulus of the two composites.

13.10 Young’s moduli of the composite systems for CNT lengths of 500 nm.

The transverse Young’s moduli of the aligned composite systems with a range of volume fractions and a CNT length of 500 nm are shown in Fig. 13.11. In contrast to the Young’s moduli in Fig. 13.10, the transverse Young’s moduli improved when the CNTs were functionalized. The enhancement is evident for CNT volume fractions greater than 10%.

13.11 Transverse Young’s moduli of composite systems for CNT lengths of 500 nm.

The primary implication of these results is that although chemical functionalization of single-walled carbon CNTs has been considered as a means to increase load-transfer efficiency in CNT–polymer composites, this functionalization has, in fact, degraded most of the macroscopic elastic stiffness components of the composite materials considered in this study. This is possibly due to the change in chemistry at the functionalization site. The carbon atom on the CNTs that is functionalized goes from an aromatic state to an aliphatic state upon functionalization. As a result, it is possible that the stiffness of the C-C bonds in the CNT is reduced at the functionalization site, thus reducing the overall composite stiffness in the direction parallel to the CNTs. However, the increase in C-C bond stiffness perpendicular to the surface of the CNT after functionalization naturally increases the overall transverse stiffness of the composite. It is important to note that these results are highly dependent on the specific materials that are considered. Different combinations of reinforcement and matrix materials would likely yield different results.

13.6 Conclusion and future trends

This chapter has discussed some of the very fundamental aspects of multiscale modeling of CNT composite materials. Although many different approaches have been taken in the literature to predict bulk properties as a function of molecular/micro-structure, all multiscale modeling approaches should follow some basic guidelines to ensure consistency with the physical assumptions of the continuum-based and non-continuum-based modeling tools. This consistency will lead to more reliable predictions and, ultimately, more efficient material development efforts.

Multiscale modeling approaches that incorporate both molecular- and continuum-level simulations are relatively new. This is mostly due to the rapid rise in computational power in the past two decades that has enabled more efficient simulations, and due to the large-scale world-wide investment in nanoscale research that was initiated around the year 2000. As a consequence, these multiscale modeling efforts are still in their infancy. There are still many aspects of multiscale modeling that will likely experience advancement in the coming years.

Many of the modeling approaches referenced herein have very limited experimental validation. This is partially due to the complicated nature of testing small-scale materials and structures, and due to the relative difficulty of fabricating test specimens of CNT composites. As multiscale modeling approaches become more sophisticated, powerful, and efficient, more reliance will be placed on their predictions for the design of engineering structures. Therefore, it is important that a greater degree of confidence must be established with the models. This will need to be achieved with more rigorous experimental validation attempts. Properly validated and benchmarked multiscale modeling tools will greatly facilitate the overall material design and development process.

Another important future direction for multiscale modeling of CNT composites is the direct modeling of time-dependent behavior. While the continuum-based tools work well for viscoelastic polymer composites whose relaxation times are on the order of seconds, minutes, hours, or days; molecular dynamics simulations of polymers typically can efficiently simulate the molecular behavior only over a span of a nanoseconds, or less. This discrepancy between time scales between tools in multiscale modeling is a fundamental road block that must be eventually cleared to truly model the time-dependent response of CNT–polymer composites.

13.7 Sources of further information

Although there have been countless journal articles in the past decade that describe multiscale modeling of nanostructured materials, the vast majority of them have focused on crystalline and other highly-ordered materials, such as metals and stand-alone carbon CNTs. By comparison, amorphous materials (e.g. polymers) are much more difficult to model and have thus received significantly less attention by multiscale modelers. Because the molecular structure of amorphous polymers lacks symmetry, simulations are restricted to much smaller systems, and the atomic potentials (force fields) are not generally as accurate as those used for crystalline systems. As a result, there are few publications outside of journal articles that address the multiscale modeling of CNT composites. Those that do address the multiscale modeling of polymer composites only do so briefly. Researchers who are interested in learning more about this field are encouraged to consult some of the journal articles referenced herein as a starting point (Frankland et al., 2003; Odegard et al., 2003; Odegard et al., 2004; Odegard et al., 2005b; Shi et al., 2005; Jiang et al., 2006; Grujicic et al., 2007; Karakasidis and Charitidis, 2007; Mokashi et al., 2007; Liu et al., 2008; Lu et al., 2008; Tserpes et al., 2008; Zeng et al., 2008; Shen, 2009; Takeda et al., 2009).

13.8 References

Aifantis, E.C. On the microstructural origin of certain inelastic models. Journal of Engineering Materials and Technology-Transactions of the ASME. 1984; 106:326–330.

Aleman, C., Munoz-Guerra, S. On the mechanical properties of poly(ethylene terephthalate): force-field parametrization and conformational analysis for the prediction of the crystal moduli. Journal of Polymer Science Part B: Polymer Physics. 1996; 34:963–973.

Bathe, K.J. Finite Element Procedures. Englewood Cliffs, NJ: Prentice-Hall, Inc; 1996.

Brown, D., Mele, P., Marceau, S., Alberola, N.D. A molecular dynamics study of a model nanoparticle embedded in a polymer matrix. Macromolecules. 2003; 36:1395–1406.

Eringen, A.C. Microcontinuum Field Theories. New York: Springer-Verlag; 1999.

Eringen, A.C. Nonlocal Continuum Field Theories. New York: Springer-Verlag; 2002.

Fan, C.F., Cagin, T., Chen, Z.M., Smith, K.A. Molecular modeling of polycarbonate. 1. force field, static structure, and mechanical properties. Macromolecules. 1994; 27:2383–2391.

Fan, C.F., Hsu, S.L. Application of the molecular simulation technique to characterize the structure and properties of an aromatic polysulfone system. 2. Mechanical and thermal properties. Macromolecules. 1992; 25:266–270.

Frankland, S.J.V., Caglar, A., Brenner, D.W., Griebel, M. Molecular simulation of the influence of chemical cross-links on the shear strength of carbon nanotube-polymer interfaces. Journal of Physical Chemistry B. 2002; 106:3046–3048.

Frankland, S.J.V., Harik, V.M., Odegard, G.M., Brenner, D.W., Gates, T.S. The stressstrain behavior of polymer-nanotube composites from molecular dynamics simulation. Composites Science and Technology. 2003; 63:1655–1661.

Griebel, M., Hamaekers, J. Molecular dynamics simulations of the elastic moduli of polymer-carbon nanotube composites. Computer Methods in Applied Mechanics and Engineering. 2004; 193:1773–1788.

Grujicic, M., Angstadt, D.C., Sun, Y.P., Koudela, K.L. Micro-mechanics based derivation of the materials constitutive relations for carbon-nanotube reinforced poly-vinyl-ester-epoxy based composites. Journal of Materials Science. 2007; 42:4609–4623.

Gusev, A.A., Zehnder, M.M., Suter, U.W. Fluctuation formula for elastic constants. Physical Review B. 1996; 54:1–4.

Herakovich, C.T. Mechanics of Fibrous Composites. New York: John Wiley & Sons, Inc; 1998.

Hutnik, M., Argon, A.S., Suter, U.W. Simulation of elastic and plastic response in the glassy polycarbonate of 4,4’-isopropylidenediphenol. Macromolecules. 1993; 26:1097–1108.

Jang, S.S., Jo, W.H. Analysis of the mechanical behavior of amorphous atactic poly(oxypropylene) by atomistic modeling. Macromolecular Theory and Simulations. 1999; 8:1–9.

Jang, S.S., Jo, W.H. Analysis of the mechanical behavior of poly(trimethylene terephthalate) in an amorphous state under uniaxial extension-compression condition through atomistic modeling. Journal of Chemical Physics. 1999; 110:7524–7532.

Jiang, L.Y., Huang, Y., Jiang, H., Ravichandran, G., Gao, H., Hwang, K.C., Liu, B. A cohesive law for carbon nanotube/polymer interfaces based on the Van der Waals force. Journal of the Mechanics and Physics of Solids. 2006; 54:2436–2452.

Jiang, M., Alzebdeh, K., Jasiuk, I., Ostoja-Starzewski, M. Scale and boundary conditions effects in elastic properties of random composites. Acta Mechanical. 2001; 148:63–78.

Jiang, M., Jasiuk, I., Ostoja-Starzewski, M. Apparent elastic and elastoplastic behavior of periodic composites. International Journal of Solids and Structures. 2002; 39:199–212.

Jin, Y., Yuan, F.G. Simulation of elastic properties of single-walled carbon nanotubes. Composites Science and Technology. 2003; 63:1507–1515.

Kang, J.W., Choi, K., Jo, W.H., Hsu, S.L. Structure-property relationships of polyimides: a molecular simulation approach. Polymer. 1998; 39:7079–17078.

Karakasidis, T.E., Charitidis, C.A. Multiscale Modeling in Nanomaterials Science. Oxford: Elsevier Science Bv; 2007.

Leach, A.R. Molecular Modelling: Principles and Applications. New York: Prentice Hall; 2001.

Li, C.Y., Chou, T.W. Multiscale modeling of compressive behavior of carbon nanotube/polymer composites. Composites Science and Technology. 2006; 66:2409–2414.

Liu, Y.J., Nishimura, N., Qian, D., Adachi, N., Otani, Y., Mokashi, V. A boundary element method for the analysis of CNT/polymer composites with a cohesive interface model based on molecular dynamics. Engineering Analysis with Boundary Elements. 2008; 32:299–308.

Lu, W.B., Wu, J., Song, J., Hwang, K.C., Jiang, L.Y., Huang, Y. A cohesive law for interfaces between multi-wall carbon nanotubes and polymers due to the Van der Waals interactions. Computer Methods in Applied Mechanics and Engineering. 2008; 197:3261–3267.

Meyers, M.T., Rickman, J.M., Delph, T.J. The calculation of elastic constants from displacement fluctuations. Journal of Applied Physics. 2005; 98:1–3.

Miranda, C.R., Tretiakov, K.V., Scandolo, S. A computational study of elastic properties of disordered systems with voids. Journal of Non-Crystalline Solids. 2006; 352:4283–4286.

Mokashi, V.V., Qian, D., Liu, Y.J. A study on the tensile response and fracture in carbon nanotube-based composites using molecular mechanics. Composites Science and Technology. 2007; 67:530–540.

Odegard, G.M. Multiscale modeling of nanocomposite materials. In: Farahmand B., ed. Virtual Testing and Predictive Modeling: For Fatigue and Fracture Mechanics Allowables. New York: Springer, 2009.

Odegard, G.M., Clancy, T.C., Gates, T.S. Modeling of the mechanical properties of nanoparticle/polymer composites. Polymer. 2005; 46:533–562.

Odegard, G.M., Frankland, S.J.V., Gates, T.S. Effect of nanotube functionalization on the elastic properties of polyethylene nanotube composites. AIAA Journal. 2005; 43:1828–1835.

Odegard, G.M., Gates, T.S., Nicholson, L.M., Wise, K.E. Equivalent-continuum modeling of nano-structured materials. Composites Science and Technology. 2002; 62:1869–1880.

Odegard, G.M., Gates, T.S., Wise, K.E., Park, C., Siochi, E. Constitutive modeling of nanotube-reinforced polymer composites. Composites Science and Technology. 2003; 63:1671–1687.

Odegard, G.M., Pipes, R.B., Hubert, P. Comparison of two models of SWCN polymer composites. Composites Science and Technology. 2004; 64:1011–1020.

Papakonstantopoulos, G.J., Doxastakis, M., Nealey, P.F., Barrat, J.L., De Pablo, J.J. Calculation of local mechanical properties of filled polymers. Physical Review E. 2007; 75:1–13.

Parrinello, M., Rahman, A. Strain fluctuations and elastic constants. Journal of Chemical Physics. 1982; 76:2662–2666.

Qi, D., Hinkley, J., He, G. Molecular dynamics simulation of thermal and mechanical properties of polyimide-carbon nanotube composites. Modeling and Simulation in Materials Science and Engineering. 2005; 13:493–507.

Raaska, T., Niemela, S., Sundholm, F. Atom-based modeling of elastic constants in amorphous polystyrene. Macromolecules. 1994; 27:5751–5757.

Raffaini, G., Elli, S., Ganazzoli, F. Computer simulation of bulk mechanical properties and surface hydration of biomaterials. Journal of Biomedical Materials Research A. 2006; 77:618–626.

Ray, J.R., Rahman, A. Statistical ensembles and molecular dynamics studies of anisotropic solids. Journal of Chemical Physics. 1984; 80:4423–4428.

Saether, E., Frankland, S.J.V., Pipes, R.B. Transverse mechanical properties of single-walled carbon nanotube crystals. Part I: Determination of elastic moduli. Composites Science and Technology. 2003; 63:1543–1550.

Shen, H.S. Nonlinear bending of functionally graded carbon nanotube-reinforced composite plates in thermal environments. Composite Structures. 2009; 91:9–19.

Sheng, N., Boyce, M.C., Parks, D.M., Rutledge, G.C., Abes, J.I., Cohen, R.E. Multiscale micromechanical modeling of polymer/clay nanocomposites and the effective clay particle. Polymer. 2004; 45:487–506.

Shi, D.L., Feng, X.Q., Jiang, H.Q., Huang, Y.Y., Hwang, K.C. Multiscale analysis of fracture of carbon nanotubes embedded in composites. International Journal of Fracture. 2005; 134:369–386.

Takeda, T., Shindo, Y., Narita, F., Mito, Y. Tensile characterization of carbon nanotube-reinforced polymer composites at cryogenic temperatures: experiments and multiscale simulations. Materials Transactions. 2009; 50:436–445.

Theodorou, D.N., Suter, U.W. Atomistic modeling of mechanical properties of polymeric glasses. Macromolecules. 1986; 19:139–154.

Truesdell, C., Noll, W. The Non-Linear Field Theories of Mechanics. New York: Springer-Verlag; 2004.

Truesdell, C.A., Toupin, R.A. The classical field theories. In: Flugge S., ed. Encyclopedia of Physics, Volume III/1: Principals of Classical Mechanics and Field Theory. Berlin: Springer-Verlag, 1960.

Tserpes, K.I., Papanikos, P., Labeas, G., Pantelakis, S.G. Multi-scale modeling of tensile behavior of carbon nanotube-reinforced composites. Theoretical and Applied Fracture Mechanics. 2008; 49:51–60.

Valavala, P.K., Clancy, T.C., Odegard, G.M., Gates, T.S. Nonlinear multiscale modeling of polymer materials. International Journal of Solids and Structures. 2007; 44:1161–1179.

Van Workum, K., de Pablo, J.J. Computer simulation of the mechanical properties of amorphous polymer nanostructures. Nano. Letters. 2003; 3:1405–1410.

Wang, Y., Sun, C., Sun, X., Hinkley, J., Odegard, G.M., Gates, T.S. 2-D nano-scale finite element analysis of a polymer field. Composites Science and Technology. 2003; 63:1581–1590.

Wang, Y., Zhang, C., Zhou, E., Sun, C., Hinkley, J., Gates, T.S., Su, J. Atomistic finite elements applicable to solid polymers. Computational Materials Science. 2006; 36:292–302.

Wu, C., Xu, W. Atomistic molecular modelling of crosslinked epoxy resin. Polymer. 2006; 47:6004–6009.

Yang, J.S., Jo, W.H. Analysis of the elastic deformation of semicrystalline poly(trimethylene terephthalate) by the atomistic-continuum model. Journal of Chemical Physics. 2001; 114:8159–8164.

Yoshimoto, K., Papakonstantopoulos, G.J., Lutsko, J.F., De Pablo, J.J. Statistical calculation of elastic moduli for atomistic models. Physical Review B. 2005; 71:1–6.

Zeng, Q.H., Yu, A.B., Lu, G.Q. Multiscale modeling and simulation of polymer nanocomposites. Progress in Polymer Science. 2008; 33:191–269.

Zhou, Z. Fluctuations and thermodynamics properties of the constant shear strain ensemble. Journal of Chemical Physics. 2001; 114:8769–8774.