## Abstract

This paper reviews and enhances numerical models for determining thermal, elastic and electrical properties of carbon nanotube-reinforced polymer composites. For the determination of the effective stress–strain curve and thermal conductivity of the composite material, finite-element analysis (FEA), in conjunction with the embedded fibre method (EFM), is used. Variable nanotube geometry, alignment and waviness are taken into account. First, a random morphology of a user-defined volume fraction of nanotubes is generated, and their properties are incorporated into the polymer matrix using the EFM. Next, incremental and iterative FEA approaches are used for the determination of the nonlinear properties of the nanocomposite. For the determination of the electrical properties, a spanning network identification algorithm is used. First, a realistic nanotube morphology is generated from input parameters defined by the user. The spanning network algorithm then determines the connectivity between nanotubes in a representative volume element. Then, interconnected nanotube networks are converted to equivalent resistor circuits. Finally, Kirchhoff's current law is used in conjunction with FEA to solve for the voltages and currents in the system and thus calculate the effective electrical conductivity of the nanocomposite. The model accounts for electrical transport mechanisms such as electron hopping and simultaneously calculates percolation probability, identifies the backbone and determines the effective conductivity. Monte Carlo analysis of 500 random microstructures is performed to capture the stochastic nature of the fibre generation and to derive statistically reliable results. The models are validated by comparison with various experimental datasets reported in the recent literature.

## 1. Introduction

Carbon nanotubes (CNTs) possess exceptional mechanical, thermal and electrical properties at an almost insignificant weight on the macroscale. Theoretically, single-walled carbon nanotubes (SWCNTs) exhibit elastic moduli of the order of 1 TPa and fracture strains of 10–30% [1–4]. These values indicate that SWCNTs have elastic moduli values that are three times those of carbon fibres and five times that of steel, at one-sixth of the weight [5]. In addition to their outstanding mechanical properties, CNTs exhibit exceptional thermal and electrical properties. The theoretical thermal conductivity of SWCNTs is usually reported as 6000 W m^{−1} K^{−1} at room temperature [6], which is approximately three times the thermal conductivity of diamond, and approximately 10^{4} times that of most polymers [7]. This value is strongly temperature-dependent, as reported by Grujicic *et al.* [8], ranging from values as high as 12 000 W m^{−1} K^{−1} at 100 K to about 2000 W m^{−1} K^{−1} at 300 K. CNTs are also found to have exceptional electrical conductivity, of the order of 10^{4}–10^{7} S m^{−1}, approximately 20 decades higher than that of most polymers [7,9]. More notably, nanotubes are capable of sustaining current densities above 10^{9} A cm^{−1} at high voltages and elevated temperatures without any change in resistance or physical deterioration [10]. This indicates that CNTs have electrical current capacity 1000 times that of copper wires [11] and at least two orders higher than those of typical superconductors [12]. Further, they have the potential to enhance polymer conductivities by as much as 10 orders of magnitude [13]. These properties, along with an aspect ratio of approximately 1000, make SWCNTs strong candidates for composite reinforcement.

Their mechanical properties make them quite appealing for use as reinforcing agents—either alone or as filler with other reinforcements, such as carbon or glass fibres. Their thermal properties make them desirable for thermal management applications, and their electrical properties make them ideal for countless other composite applications, such as electromagnetic shielding, electrostatic charging prevention and damage sensing [14,15]. The very low density and high conductivity of CNTs can convert an insulating polymer into an electrical conductor without adversely impacting the other desirable properties of the thermoplastic. Collectively, all of these properties make nanotube composites of strong interest to many various fields, including the aerospace, automotive, medical and military sectors.

Because of the interest in these new composites, a persistent need to predict their physical properties has developed. Several researchers have sought to characterize these materials experimentally. Skákalová *et al.* [16] measured the electrical conductivity, stress–strain characteristics and Raman spectra of SWCNT composites of various volume fractions. Zhu *et al.* [17,18], Li *et al.* [19] and Sun *et al.* [20] all measured the elastic moduli and stress–strain characteristics. In addition, Huxtable *et al.* [21], Xu *et al.* [22], Wang *et al.* [23], Moisala *et al.* [9] and Cai & Song [24] all focused on the thermal properties of SWCNT composites.

While experiments have already demonstrated the utility of SWCNTs in composite materials, the need to accurately predict properties of nanocomposites using a proper model remains critical because the sheer number of different possible composites limits the experimental practicality. Owing to the copious combinations of CNT and polymer constituent phases, numerical and analytical methods can help reduce the number of expensive and time-consuming experiments. Nanocomposites can vary by nanotube type, polymer matrix, processing method, aspect ratio, volume fraction and nanotube alignment. A review of several hundred papers on the topic of polymer nanocomposites clearly revealed a large variation in the reported values of electrical properties as a function of CNT type, polymer matrix and processing method [25]. Athreya *et al.* [26] discovered that the electrical conductivity of a nanocomposite can vary by as much as five orders of magnitude as a result of different formation methods. Further, a percolation model developed by Bao *et al.* [27,28] demonstrated that the electrical conductivity depends upon the dispersion and alignment of the nanotubes. Thus, this manifold of controllable material parameters needs to be taken into account and studied.

Because of the aforementioned needs, many researchers have sought to model these new nanocomposites. Tserpes *et al.* [29] developed a multi-scaled modelling approach to test the tensile behaviour of composites in which nanotubes are modelled independently and then approximated as a beam element, and integrated into a representative volume element (RVE) that can next be treated using finite-element analysis (FEA). Li & Chou [30] presented another multi-scale approach to test their compressive behaviour. Xiao & Gillespie [31] reported a nanomechanics model for predicting the elastic properties of SWCNTs as composite reinforcement and then used a micromechanics model to calculate the elastic properties of the composite as a whole. Shi *et al.* [32] used a micromechanics model as well, and they included waviness and agglomeration in the calculation of elastic properties. Song & Youn [33] calculated the effective thermal conductivity of nanotube/polymer composites using the control volume finite-element method (FEM). All of these models used a single CNT that traversed the entire RVE. They assumed that the effective properties of the composite as a whole can be calculated based on a tiny RVE with a single nanotube for reasons of computational simplicity.

Besides the models that calculate effective properties based upon a single nanotube RVE, several models did not make this assumption. Xue [34] presented a numerical method based on Maxwell theory to calculate the thermal conductivity of nanotube/oil and nanotube/decene composites. Bagchi & Nomura [35] approximated nanotubes as spheroidal inclusions, and used effective-medium theory to calculate the thermal conductivity of aligned nanotube/polymer composites. One of the most versatile models is presented by Zhang & Tanaka [36]. Using the hybrid boundary node method in conjunction with the fast multipole method, they calculated the thermal properties of CNT composites, while including waviness and alignment for a variety of nanotube configurations. Odegard *et al.* [37] modelled the nanotubes as an effective continuum fibre using the equivalent-continuum modelling method along with a micromechanics approach to calculate mechanical properties of CNT/polymer composites. Seidel & Lagoudas [38] modelled nanotubes using a composite cylinders approach, and used the Mori–Tanaka method [39] to obtain effective elastic properties for the composite as a whole. They considered both aligned and randomly oriented nanotubes and also consider nanotube bundling. Spanos & Kontsos [40] presented one of the very early approaches for determining elastic properties of nanocomposites using the method of Monte Carlo simulation [41].

These efforts reflect many of the great leaps made in understanding these new nanocomposites. However, to date, scarce modelling information is available for the nonlinear properties of either CNTs or the polymer matrix for a realistic composite size. Elastic nonlinearity must be considered to accurately model the stress–strain relationship of these composites. Also, the effect of the nonlinear thermal conductivity of SWCNTs on the effective thermal conductivity of the composite is of interest. The model discussed herein studies the effects of these two types of nonlinearity on the physical properties of SWCNT-reinforced polymer composites using a novel method described by Esteva [42] and Spanos & Esteva [43].

A proper model is also needed to predict the electrical characteristics of nanocomposites. Several numerical models provide good representations of nanocomposites. Pike & Seager [44] were the first to study percolation in random two-dimensional stick models using Monte Carlo simulations. One of the original resistor network models was developed by Kirkpatrick [45] as an extension to established lattice models. Both two- and three-dimensional disordered stick systems were simulated by the Monte Carlo method. More recently, Behnam & Ural [46] used Kirchhoff's laws and Monte Carlo simulations to calculate the electrical properties of straight, randomly oriented CNTs in stacked two-dimensional layers. Dalmas *et al.* [47] presented a three-dimensional model of wavy, randomly oriented nanotubes using the resistor network method. They included constant contact resistance and solved the model using the FEM. Hakobyan *et al.* [48] also used the FEM to determine the effective conductivity and percolation threshold of a two-dimensional lattice of highly conductive triangles. Foygel *et al.* [49] studied only the percolation threshold with their spanning network model of straight, randomly oriented and soft-core CNTs. Monte Carlo, statistical percolation and resistor network methods were used by Hu *et al.* [50] to determine electrical conductivity and thresholds in polymer nanocomposites. The pseudo-three-dimensional model of Jack *et al.* [51] was based on Kirchhoff's current law (KCL) and Monte Carlo simulations; it included stochastic nanotube parameter inputs, properly scaled sample dimensions and periodic geometry on the boundaries. Li & Chou [52] presented a method of representing contact resistances and used their direct electrifying algorithm to calculate effective conductivity and identify the backbone [53]. The three-dimensional model of Sun [54] used the resistor network method and calculated contact resistances to find the effective conductivity and percolation threshold. Topinka et al. [55] also used the resistor network method combined with KCL and constant contact resistances to determine the effective conductivity and visualize the currents and voltages in the CNT networks.

Each of these models provides important insights into the electrical conduction behaviour of nanocomposites. The vast majority used a combination of Monte Carlo and resistor network methods. Pike & Seager [44] were the first to establish Monte Carlo simulations as the primary method of analysing random percolating systems. Kirkpatrick [45] was the first to prove the viability of a resistor network model in determining the electrical properties of a system of random sticks. Though they used similar methods, all the models were unique in how details such as morphology, contact resistance, CNT morphology and CNT resistance were represented. Only a few models used periodic geometric boundaries. By relocating the parts of fibres extending outside the sample to the opposite boundary, continuity between sample structures was maintained, insulating boundary conditions were validated, and the structure was representative of an infinite system [45,51]. Also, contact resistance between CNTs is one of the most critical physical aspects of the problem and must be incorporated into the model [54]. Some numerical models involved broad assumptions about contact resistance or ignore the tunnelling effect altogether. Several authors [14,45,55,56] emphasized the impact of the number of CNT contact points on the results. For this reason, it is important for models to incorporate nanotube waviness and realistically represent the number of contact points. The other major factor to include in a numerical model is the intrinsic nanotube resistance [14]. Of the models detailed above, none incorporated all four of these important factors. Therefore, a new model is necessary to simultaneously predict percolation, identify conducting backbones and estimate the electrical conductivity. These are indeed the principal motivating factors for developing a comprehensive and versatile numerical method of accurately representing the electrical behaviour of CNT-based composites.

## 2. Elastic and thermal properties

### (a) Methodology

The discussed nonlinear FEM model for thermal and elastic properties was first reported in Elsbernd & Spanos [57]. It is based upon the embedded fibre FEM (EFFEM), which was first reported in Esteva [42] and Spanos & Esteva [43]. In this method, an RVE of a user-specified volume fraction of SWCNTs is developed in three distinct steps. First, the nanotube morphology throughout the RVE is generated in a random manner, incorporating nanotube waviness, entanglement and CNT diameter/length distributions. The RVE thickness is assumed to be much smaller than the width or diameter, and thus the morphology is only generated in two dimensions. By making this assumption, the problem becomes a two-dimensional plane-stress problem, greatly reducing the computational effort required. Once the morphology is generated, a square finite-element mesh size is chosen, and the nanotubes are divided into two new shorter nanotubes when they cross the boundary from one element into another. This process is known as partitioning and is necessary to perform the final step of the method—fibre embedment. After partitioning, all nanotubes lie completely inside elements. This attribute allows the direct addition of their elastic or thermal properties to those of the parent element in a finite-element stiffness matrix sense. This direct combination of properties is known as fibre embedment. These three steps form the foundation upon which a nonlinear finite-element solution approach can be used to determine a nonlinear stress–strain curve and thermal conductivity for any arbitrary SWCNT/polymer composite.

#### (i) Representative volume element generation

Unlike the majority of computational models that currently exist for SWCNT/polymer composites, the EFFEM does not generate the composite for analysis based upon smaller RVEs with a single nanotube that traverses the entire RVE. Rather, by directly incorporating nanotube properties into the matrix, a complex nanotube morphology where some nanotubes overlap and curve can be modelled. With this in mind, a scheme for generating complex nanotube morphologies is developed. In this scheme, each nanotube is generated individually, and the length and diameter of each nanotube are randomly generated from Weibull and lognormal distributions, respectively. The details of this formulation can be found in Esteva [42] and Spanos & Esteva [43]. Using this method, complex RVE microstructures can be generated. One particular such microstructure is shown in figure 1.

#### (ii) Fibre partitioning

After generating the microstructure, the next step is to partition the nanotubes that cross from one element into another. This step is necessary because, in order to combine the stiffness matrices of the fibres within an element with the stiffness matrix of the polymer matrix for that element, the reinforcing fibre must lie completely within the element. To accomplish this partitioning, the coordinate pairs for every nanotube segment are examined to see whether they lie in the same element. If they do not, a line equation is generated for the segment. This equation is compared with the equations for the lines of the RVE mesh grid and the points of intersection are found. The nanotube is then divided at these points. For a more detailed description of the partitioning process, again see Esteva [42] and Spanos & Esteva [43]. A section of a partitioned microstructure section is presented in figure 2.

#### (iii) Embedded fibre method

Once the nanotube morphology has been generated and the fibres are partitioned, the final step to prepare the RVE for traditional FEA is to ‘embed’ the nanotubes in the polymer matrix. The process of embedding, as described above, is to directly add the nanotube ‘stiffness’ properties to those of the surrounding matrix. This is accomplished for each element in the RVE finite-element mesh. Regardless of the type of problem to be solved (elastic, thermal, electrical, etc.), the same embedding process is applied to the RVE.

In general, there are three steps to the embedding process. The first step is to use the traditional finite-element formulation to derive the stiffness matrix of the polymer element. This formulation can be found in any basic finite-element text, such as those by Akin [58], Desai & Abel [59] and Zienkiewicz *et al.* [60].

The second step is to calculate the stiffness matrix of each nanotube within the element. In the model discussed herein, the nanotubes are approximated as line elements. This approximation is accurate and computationally efficient, as shown by Konrad & Graovac [61].

The final step in the embedding process is to add the nanotube stiffness properties to those of the surrounding polymer for the element. This process is simple, accomplished by direct addition, and is encapsulated by the equation
2.1where is the equivalent stiffness matrix, *K*^{e}_{elem} is the stiffness matrix of the polymer in the element, and is the sum of the stiffness matrices of all the nanotubes in the element.

By approaching the problem in this way, the need for complicated, manual meshing around the nanotube–polymer matrix is eliminated. A simple, square finite-element mesh is all that is necessary for the EFFEM, and this automation of the model, without human intervention, directly follows from the process. This is the primary advantage of this approach. However, as in most models, an important assumption is made. The EFFEM assumes perfect bonding between the nanotubes and the polymer matrix. While initial experiments indicated that this may not be a valid assumption, recent findings have found this assumption justified in two ways. First, dramatic advancements have been achieved in interfacial bonding through chemical functionalization of the nanotube sidewalls and endcaps. Garg & Sinnott [62] reported that even a high degree of sidewall functionalization will degrade the mechanical strength of SWCNTs by only 15 per cent. Demonstrating the effectiveness of functionalization, Zhu *et al.* [17] reported that, using acid treatment and fluorination of the nanotubes, they are able to achieve a 30 per cent increase in elastic modulus of an epoxy composite with the addition of only 1% weight fraction (1 wt%) of SWCNTs. Second, conventional observation of reinforcing inclusions indicates that interfacial weakening often occurs, causing a detrimental effect on the load transfer properties of the composite. Esteva & Spanos [63] showed that, for small volume fractions, less than or equal to 30 per cent, the interface weakening effect of elliptical inclusions was not relevant for nanotube-reinforced composites. Based on these two justifications, the EFFEM is a valid approach to modelling nanocomposites and this fact is supported by the results of this paper.

### (b) Nonlinear model

After establishing the EFFEM as the framework in which to approach the problem of modelling an SWCNT-reinforced polymer composite while incorporating many realistic and essential features of such a material, the next step is to determine how to model nonlinearity. The nonlinear mechanical properties of SWCNT-reinforced composites for small strains are well observed in several experiments. Sun *et al.* [20], Zhu *et al.* [17], Skákalová *et al.* [16] and Gojny *et al.* [64] all reported clearly nonlinear stress–strain relationships for nanotube-reinforced composites at various volume fractions. Because strain levels are small, the nonlinearity was dominated by material properties rather than geometrical considerations. The nonlinear thermal properties of composites is a less explored area. Xu *et al.* [22], Hong & Tai [65], Cai & Song [24] and Biercuk *et al.* [66] investigated this area with mixed results. However, for this nonlinear analysis, the deformation of the finite-element mesh is not an issue; thus, only material nonlinearity is important.

Because material nonlinearity is the primary concern, two different approaches are considered to incorporate the nonlinear properties of both the CNTs and the polymer matrix. The first is the incremental approach. By solving a linear problem in a series of small increments, the nonlinearity of the materials can be incorporated in the model through a recursive updating procedure. The second approach is an iterative technique. In this approach, the full linear problem is solved repeatedly, again updating the properties of the nanotubes or the polymer after each step. The incremental approach is more versatile, yet more computationally expensive. It allows the user to recover intermediate solution values, which are useful for calculations such as a stress–strain curve. The iterative approach is more robust and usually less computationally costly. However, the iterative approach provides only final solution values; therefore, it has limited applicability for some applications. The general concepts of both approaches are described in detail by Desai & Abel [59]. To decide which approach to use for the elastic analysis and which to use for the thermal analysis, an understanding of the mechanisms of nonlinearity in both CNTs and epoxies is needed.

#### (i) Nonlinear properties of carbon nanotubes

The nonlinear elastic and thermal properties of CNTs are well documented in the current literature, both theoretically and experimentally. For the nonlinear elastic properties, Tserpes *et al.* [29] presented a detailed finite-element model of single nanotubes using a structural mechanics approach. They identified the stress–strain relationship for armchair and zigzag SWCNTs. Additionally, Tserpes *et al.* [3] presented a progressive fracture model to calculate the stress–strain relationship for SWCNTs with perfect structures and for those with imperfections. Belytschko *et al.* [4] developed a molecular mechanics model to simulate nanotube fracture. They reported the stress–strain curve for perfect and imperfect zigzag nanotubes as well as a curve for various different perfect nanotube structures. Natsuki *et al.* [67] described a structural mechanics approach to modelling nanotube elastic properties. Finally, Meo & Rossi [2] presented a molecular-mechanics-based FEM for the prediction of SWCNT elastic properties. Each of these models used a slightly different approach, but they all identified approximately 10–15% strain as the region where significant nonlinear behaviour is first exhibited by SWCNTs.

The nonlinear thermal properties of SWCNTs have also been modelled, though less frequently than the elastic properties. Grujicic *et al.* [8] presented a molecular-dynamics-based model of the thermal conductivity of SWCNTs in which they investigated the nonlinear dependence of thermal conductivity on temperature. The strong dependence they reported provides, partially, the motivation for the nonlinear thermal analysis performed in this paper.

Kawamura *et al.* [68] also reported the temperature dependence of SWCNT thermal conductivity using a molecular dynamics simulation. Further, Xu *et al.* [22] reported this effect indirectly by providing the temperature dependence of the thermal conductivity of composites reinforced with SWCNTs.

Along with the theoretical modelling evidence presented above, experimental evidence corroborates the nonlinear thermal and elastic properties of SWCNTS. Yu *et al.* [69,70] performed tensile load testing using an atomic force microscopy tip to document the nonlinear stress–strain relationship of SWCNTs. Their results showed the onset of nonlinear behaviour to occur at greater than 10 per cent strain. Walters *et al.* [71] also reported the lower bound of the yield strain of SWCNTs to be approximately 6 per cent. For thermal properties, Shaikh *et al.* [72] reported a range of thermal conductivities for CNT ropes to range from 0.1 to 2000 W m^{−1} K^{−1}.

#### (ii) Nonlinear properties of epoxies

Epoxy composites exhibit nonlinear elastic properties at a much lower strain level. Zhu *et al.* [17,18] presented stress–strain curves for neat epoxy in which nonlinearity was shown at 3 per cent strain. Li *et al.* [19] presented similar effects with a load–displacement curve of neat epoxy. Sun *et al.* [20] and Gojny *et al.* [64] also indicated 3 per cent strain as the onset of nonlinearity in neat epoxy. Each group indicated epoxy failure at approximately 6 per cent without reinforcement. The significantly lower onset of nonlinear behaviour and failure in the epoxy than in the lone nanotubes led to the conclusion that the nonlinear behaviour of the nanotubes may be neglected and that of the epoxy must be taken into account for an accurate elastic model at small strains less than 6 per cent.

In general, polymers were reported to have constant values of thermal conductivity. Xu *et al.* [22] provided a constant thermal conductivity of approximately 0.2 W m^{−1} K^{−1} for their poly (vinylidene fluoride) (PVDF) polymer. Wang *et al.* [23] reported a constant thermal conductivity of approximately 0.18 W m^{−1}K^{−1} for the neat epoxy resin used in their experiments. Moisala *et al.* [9] reported a constant value of approximately 0.255 W m^{−1} K^{−1} for epoxy consisting of a bisphenol-A resin and an aromatic amine hardener. Du *et al.* [7] also reported a constant thermal conductivity of approximately 0.18 W m^{−1} K^{−1} for their poly(methyl methacrylate) (PMMA) polymer. Finally, Cai & Song [24] used a semicrystalline polyurethane (PU) dispersion with a constant thermal conductivity of approximately 0.15 W m^{−1} K^{−1}. Because of the lack of any indication of a nonlinear temperature dependence of the thermal conductivity of the neat epoxy, only the nonlinear thermal conductivity of the nanotubes should be considered for a thermal model of CNT-reinforced polymer composites.

In the present paper, the nonlinear elastic properties of the epoxy matrix are considered in the elastic analysis of the model, and the nonlinear thermal properties of the CNTs are considered for the thermal analysis. These choices are based upon the understanding of the nonlinear thermal and elastic properties of both CNTs and polymer matrix materials. With these choices made, the next step is to select the nonlinear approach to take in each case.

#### (iii) Choice and implementation of nonlinear approaches

After careful consideration of the incremental and iterative methods, the incremental technique is chosen for the elastic analysis, and both techniques are used for the thermal analysis. The incremental technique is chosen both for its versatility and because it most closely resembles the physical phenomenon of a tensile strain test. The iterative technique is not applicable in the elastic analysis of the current model because it does not provide any information about the stresses and strains at intermediate load increments. Thus, a stress–strain relationship cannot be obtained from an iterative approach in the elastic case. In addition, the availability of experimental data for comparison contributed to this choice.

For the thermal analysis, both an incremental approach and an iterative approach are adopted for comparison. An iterative approach is more readily used than an incremental approach and is often more quickly convergent than the incremental approach. For this reason, it is used to calculate an effective thermal conductivity for the composite. However, it does not provide any conductivity information at intermediate temperature increments. Thus, an incremental approach is also used.

After deciding which approach to take for the thermal and elastic analyses, a nonlinear finite-element model is developed, using the embedded fibre method as a foundation. After generating the nanotube geometry, partitioning the fibres and embedding the fibres, the stiffness matrix of the entire RVE is assembled using the standard finite-element procedures. The standard force balance equation for FEA is
2.2where **K** is the system stiffness matrix assembled from the element stiffness matrices defined in equation (2.1), ** u** is the displacement or temperature vector for all of the mesh nodes in the RVE and

**is the reaction force or temperature gradient vector at the nodes. For the elastic analysis case, the general finite-element problem is solved using a displacement-controlled approach, rather than force-controlled, because a uniform displacement on an entire edge is needed to generate a stress–strain curve for the RVE as a whole. To use the incremental solution approach, displacement increments are applied to the elastic problem, and temperature increments are applied to the thermal problem. After each increment is applied, the resulting displacements or temperatures of every node in the mesh are calculated.**

*F*For the elastic case, there are two degrees of freedom (d.f.) per node, i.e. the *x*- and *y*-displacements. With the displacements at each node in hand after an incremental displacement is applied, the strain at each node can be calculated and interpolated to the Gaussian quadrature points using the strain matrix **B** used in FEA texts and the equation
2.3

In two dimensions, the strain vector ** ε** is a three-component vector including the strain in the

*x*- and

*y*-directions,

*ε*

_{x},

*ε*

_{y}, as well as the shear stress,

*γ*

_{xy}. With the strain, the stress at each Gauss point can be calculated using the material matrix

**D**and the relation 2.4

After finding the strain vector at a node using equation (2.3) and using it to find its corresponding stress vector by equation (2.4), the von Mises stress at each Gauss point in each element can be calculated. The von Mises stress is then used to refer to the stress–strain curve of the neat epoxy to calculate the tangent modulus at any particular Gauss point. In this manner, the modulus at each Gauss point throughout the RVE can be updated after each displacement increment. This procedure allows the nonlinearity of the material to be effectively incorporated into the model.

For the thermal problem, the incremental approach follows essentially the same approach as for the elastic problem, with only a few differences. The first difference between the thermal and elastic problems is that the thermal case requires only 1 d.f. per node, i.e. the temperature. This in turn greatly reduces the computational cost to solve the problem because the number of equations to solve is halved. The second difference is that, instead of updating the elastic modulus of the material, the thermal conductivity of the nanotubes is updated after each incremental temperature increase. This process is simply performed with the knowledge of the temperature at each of the nodes throughout the RVE. Using the shape functions of a quadrilateral element, the temperature is interpolated from the nodal temperatures to define a temperature at each of the nanotube endpoints. Because the nanotube thermal conductivity is quite high—approximately four orders of magnitude greater than that of the polymer matrix—the average of the two nanotube endpoint temperatures is assumed as the constant temperature along that nanotube. The thermal conductivity of that particular nanotube is then updated based upon the curve of Grujicic *et al.* [8]. The equation used to update the thermal conductivity is given by
2.5where *K*_{nt} is the thermal conductivity of the nanotube and *T* is the temperature of the nanotube in degrees celsius. This procedure is repeated for every nanotube in each element and after each incremental temperature step. In addition to the incremental approach, the iterative technique is also used for the thermal problem. In this process, the finite-element solution is solved with a single, large temperature difference. Using the temperature data recovered from this solution, the thermal conductivities are updated using the same procedure as the incremental approach. The problem is then solved again using these data. This process is repeated until the solution converges to a solution.

### (c) Results and discussion

To calculate an effective stress–strain curve of the composite, a tensile displacement is applied for the elastic problem and transverse movement of the structure is unrestricted. By adopting these boundary conditions, transverse stresses at the bottom nodes are eliminated and a direct correlation can be made between the applied displacement and the stresses recovered at the bottom nodes that are held in the tensile direction.

For the thermal problem, the bottom edge of the RVE is held at 0^{°}C and the top edge is held at 100^{°}C. The sides of the RVE are insulated and the effective thermal conductivity of the RVE is calculated using the sum of the heat fluxes through the bottom nodes held at 0^{°}C.

To qualitatively verify the EFFEM method, a preliminary elastic and thermal analysis is performed. To verify the elastic results, the 1 μm^{2} RVE is displaced to a strain of 3.5 per cent for an arbitrary nanotube geometry of 1 wt% SWCNTs. By visualizing the displacement contours in the tensile stretching direction along with the nanotube geometry, the reinforcing effect of the nanotubes can clearly be seen by noticing the areas of the plot in which colour intensity remains constant. The true power of the reinforcement is seen by examining the displacements at a high strain value, after the polymer has begun to exhibit a large degradation of elastic modulus. These two contours are shown in figure 3 for comparison.

For the thermal problem, the boundary conditions above are used and a thermal contour plot is produced. The temperature contour for one particular nanotube geometry of 1 wt% SWCNTs is presented in figure 4. Again, the enhancement of thermal conductivity can be seen by noticing the areas where the colour intensity remains constant.

After confirming the general validity of the EFFEM for both the elastic and thermal nonlinear problems, the model is compared with several experiments by using the same model parameters that are reported in the corresponding experimental reports.

#### (i) Elastic results

Before proceeding to comparison with experimental data, a Monte Carlo mesh convergence analysis (MCMCA) is performed to determine the mesh refinement and incremental step size necessary for reliable results. Esteva [42] performed MCMCA for the linear EFFEM model and varied the number of random microstructures used; it was determined that 500 microstructures are sufficient to derive statistically reliable values of effective moduli of elasticity and thermal conductivities. This same number of microstructures is used in the present model. From the MCMCA, it is determined that 60 divisions and 100 incremental steps are sufficient for an accurate finite-element model that captures the nonlinear elastic properties of the polymer.

The first experiment chosen for comparison is that of Sun *et al.* [20], and the stress–strain curves are shown in figure 5. The discrepancy between the two curves is no more than 5 per cent at any point along the curve and is explained by two factors. First, an increase in mesh refinement of the model would bring the model curve closer to the experimental curve. Because the aim herein is merely to show that the model is accurate, not exact, the mesh refinement is not changed. Second, the model curve is an average over 500 random microstructures, whereas the experimental curve is a single realization of one particular nanotube geometry. Any single geometry curve strongly depends on the particular orientation of nanotubes in that sample, as is illustrated by comparing the maximum and minimum curves of the 500 microstructures to the same experimental curve in figure 6.

The second set of experimental data chosen for comparison is that of Zhu *et al.* [17]. They also used 1 wt% SWCNTs in an epoxy matrix. They performed a tensile test to 5 per cent strain and report a different stress–strain curve for their neat epoxy. A comparison of the stress–strain curves of the model and the experiment is performed, and is shown in figure 7.

Akin to the first experiment, the model is accurate to within 5 per cent for the entire curve. Again, this feature could be improved by an increase in mesh refinement. These two comparisons show that the proposed nonlinear model effectively captures the material nonlinearity of polymer-based nanotube composites.

#### (ii) Thermal results

As a parallel study, the effect of nonlinear thermal properties on the effective thermal conductivity of the composite is also examined. To model the nonlinear thermal conductivity of SWCNTs, both incremental and iterative approaches are used. In both approaches, the thermal conductivity of the composite is calculated using the heat fluxes at the bottom nodes after the final step. The equation used to calculate the effective thermal conductivity is
2.6where is the sum of the heat fluxes across the bottom nodes, *A* is the area of the RVE, *t* is the thickness of the RVE and *ΔT* is the temperature difference between the bottom and top of the RVE. To determine which nonlinear approach is more suitable, as well as how many divisions are necessary in the mesh, MCMCA is again performed. It is found that 130 divisions and 10 incremental steps are sufficient for the incremental approach. For the iterative approach, it is also found that 130 divisions are sufficient. However, only two iterative steps are needed to capture the nonlinearity of the thermal conductivity of the nanotubes. Based on these findings, the iterative approach is selected because it converges to a single value after only two steps. A mesh refinement of 130 divisions is also chosen because there is only 1 d.f. per node in the thermal analysis, and thus few iterative steps are needed.

Next, the thermal results are compared with various experimental datasets. Because many of the experimental results are reported in weight fractions, a conversion from weight fraction (*wf*) to volume fraction (*vf*) is performed using the equation [73]
2.7A matrix density of 1.2 g cm^{−3} is assumed because most studies do not detail the density of their polymer, and this value is typical for epoxies. Further, provisions are made to account for the high thermal resistance between the nanotube–matrix and nanotube–nanotube interfaces. Xue [74] found that the effects of the interfacial resistance, or the Kapitza resistance, between the nanotubes and the surrounding matrix can effectively reduce the thermal conductivity of the nanotubes by several orders of magnitude. Using the work of Huxtable *et al.* [21], Wilson *et al.* [75] and Xue [74] as justification, a phenomenological degradation of nanotube thermal conductivity by a factor of 20 is incorporated into the model.

The data of Wang *et al.* [23], Moisala *et al.* [9], Hong & Tai [65], Biercuk *et al.* [66], Cai & Song [24] and Xu *et al.* [22] are selected to provide a broad spectrum of polymers and volume fractions for comparison. A comparison with these six experiments is presented in figure 8. Equation (2.7) and the reported matrix thermal conductivity for each experiment are used in this comparison.

For every comparison, except that to Xu *et al.*, the thermal conductivity predicted by the model is within approximately 33 per cent of the value reported by the corresponding experiment. While this value may seem large at first, this agreement is good after considering the poorly understood and complex mechanisms that inhibit efficient heat conduction; the assumption of perfect bonding between the matrix and the nanotubes for all values of volume fractions is also subject to further scrutiny.

## 3. Electrical properties

### (a) Methodology

Next, an electrical properties model is adopted. First, the morphology of a CNT-based composite material is generated in the same manner as the elastic and thermal properties model. Then, a spanning network algorithm determines the existence of a complete percolating nanotube network. For simplicity, the morphology of the nanocomposite is represented by a pseudo-three-dimensional RVE. The RVE is a two-dimensional square with an assumed thickness that approximates a three-dimensional structure. The CNT fibres are randomly placed within the RVE at the desired concentration or volume fraction. Next, each CNT is given extra nodes as searching points for the spanning network algorithm. Then, this algorithm uses the specified bond criteria to determine the connections between fibres. Finally, the algorithm determines whether a complete network of CNTs spans the longitudinal direction of the RVE.

### (b) Spanning network identification

After creating a realistic nanocomposite morphology in each RVE, the microstructure must be analysed to determine the existence of a complete spanning CNT network. A complete spanning network allows an electrical connection to be made between the top and bottom edges of the sample structure. This percolation condition is the determining factor in effective electrical properties of nanocomposites. Without percolation, no electrical current flows through the RVE.

The RVE is divided into a regular grid of searching bins. Each nanotube is assigned to a bin based upon its location in the RVE. The goal of using bins is to accelerate the searching process. A visualization of a typical bin is shown in figure 9.

Since the fibres are assembled into bins, the remainder of the spanning network algorithm depends on the bonding criterion. This criterion determines whether any pair of points, each from a different CNT, is connected [44]. If multiple fibres and bonds belong to the same chain, they form a CNT network. A complete spanning network consists of an unbroken series of fibres and connections that touch the top and bottom of the RVE. Theoretically, when this happens, the RVE is percolated and electrical current may flow in the nanocomposite. The model's criterion states that a connection is established between two nodes if the position of a fibre node is within the searching region of the fibre being checked. Connections between all combinations of fibres in the RVE are included in the search. Pike & Seager [44] proved that there is little difference between using circles or squares as searching regions; squares are used herein. The searching regions in this model are variable. The size of each variable searching region (VSR) depends on the diameter of each fibre to which the pair of nodes belongs plus the electron tunnelling distance. Figure 10 illustrates how the VSR is used for fibres with different diameters.

After the fibres are sorted into bins and the bond criterion is established, extra nodes are added along the fibre segments. These nodes are used as additional searching points for determining connections between fibres. At this stage, nodes exist only at the ends of each CNT segment in the RVE. Adding additional nodes creates more VSR searching points, but also increases the number of node pairs that must be checked. Therefore, it is critical that the new nodes are properly spaced to minimize the computational cost and to simultaneously account for the effective tunnelling region around each fibre. New node spacing is determined by the inequality
3.1where *n* is the number of new nodes to add, *l*_{seg} is the original fibre segment length, *D*_{min} is the smallest possible fibre diameter, and *TunL* is the electron tunnelling distance. The VSR boxes will leave no gaps along the fibre, ensuring all connections are accounted for without increasing the computational cost of the model more than necessary.

Finally, the spanning network algorithm searches sequentially for CNT clusters in the RVE. This is accomplished in a series of steps, beginning with determining if the random RVE has a potential spanning cluster. A complete CNT network only exists if at least one fibre is in contact with the top boundary and at least one other fibre is in contact with the bottom of the RVE. After all top and bottom edge fibres are identified, the top fibres are used as the building blocks for the CNT network. The VSR criterion is used to find other fibres connected to the top fibres. The nodes of the top fibres are the searching points for establishing these new connections. This part of the model is then repeated and the newly found fibres become the searching points. With each new iteration, the fibres being added to the network approach the bottom boundary. If after a particular iteration no new fibres satisfy the bonding criterion, this part of the algorithm stops. Figure 11 illustrates how the network is built in the model.

After the iterations are complete, the connections are checked for duplicates. Each fibre can have multiple connections with another fibre, but only one connection is made between the same two segments. This means CNTs that lie parallel to each other or cross multiple times are modelled as parallel resistors. Finally, the percolation status of the RVE is determined. If any fibre in the network is identified as touching the bottom edge, then a complete network of CNTs exists in the nanocomposite. Therefore, each RVE with a percolating network has the potential to conduct electrical current along the longitudinal direction, which is determined by the resistor network portion of the model.

### (c) Equivalent resistor network

An important component of the resistor network model is the representation of fibres and contacts as electrical resistor elements. It is critical that both the fibre and contact resistances are modelled accurately because they have a significant impact on the results. Commonly, electrical resistance is calculated by the equation
3.2where *l* is the length of the CNT or CNT segment, *d* is nanotube diameter and *σ* is CNT intrinsic conductivity [50,54,76]. Equation (3.2) is adopted in the present model for its simplicity and accurate representation of resistance. The nanotube resistances calculated for the model encompass a range from 0.16 to 52.33 kΩ with a mean of 5.89 kΩ.

With all the fibres accurately represented as electrical resistors, the next step is to model the contact points in the resistor network. The contacts between CNTs throughout the RVE are determined by the spanning network algorithm. These connections are converted to equivalent resistors to accurately reflect electrical transport behaviour at CNT–CNT junctions. Several authors [54,76–80] have reported that the electrical properties of nanocomposites are governed by electron hopping or tunnelling. Experiments and models revealed that the electron tunnelling results in nonlinear current–voltage behaviour in nanocomposites [81,82]. The resistance between nanotubes is a combination of direct contact resistance and tunnelling resistance [83]. Direct contact resistance was reported to be in the range of 100–400 kΩ [84,85]. The method incorporated into the present model is taken from the work of Li *et al.* [76] and Li & Chou [52]. Li *et al.* [76] used the Simmons formula to calculate current density and hence tunnelling resistance between crossed CNTs with a finite contact area. Theoretically, the tunnelling resistance reached 10^{19} Ω at a tunnelling distance limit of 1.8 nm. This 1.8 nm limit is used in the present model as the maximum polymer distance through which an electron can tunnel from one CNT to another. The tunnelling distance can vary from 0 nm to 1.8 nm. Therefore, the distances for all connections in the model are sampled randomly from a normal distribution similar to the method of Li & Chou [52].

For each connection established in the resistor network, a random tunnel distance is generated from the assumed normal distribution. Then, this value is used to calculate a tunnelling resistance from an exponential fit to the data reported by Li *et al.* [76]. The lower limit of contact resistance is set equal to 10^{5} Ω to reflect direct contact. Maximum resistance is limited to be 10^{16} Ω to represent the highest possible tunnelling distance and to prevent singularities in the solver. Next, each CNT and connection in the spanning network is assigned a resistance. An example of connected fibres in the RVE is shown schematically in figure 12.

After the complete spanning network is converted to an equivalent resistor network, the next step is to use KCL in conjunction with FEA to solve for the system voltages and currents.

### (d) Details of the model

The application of KCL to resistor networks is well established. Kirkpatrick [45] was among the first to use KCL on two- and three-dimensional site and bond percolation models based on resistor networks. The foundational principle of KCL is conservation of electrical charge. It states that the sum of the currents at a given point, or in this case node, must equal zero. Also, the total current entering and exiting the nanocomposite system must be the same. Essentially, the movement of electrons through the percolating model can be viewed as a fluid flowing through the resistor network.

The two pieces of information required to apply KCL and FEA are the calculated resistances and nodal connection information. Each component in the network is modelled as a one-dimensional linear resistor element. This includes both the nanotubes and the tunnelling contacts. Applying Ohm's law at a node yields the equation 3.3which can be represented in matrix form as 3.4where the element coefficient matrix is 3.5

Equation (3.4) reflects the finite-element formulation for a single resistor element. The nodal currents included in this equation are with respect to the particular resistor element of interest. The complete CNT network is merely a system of linear algebraic equations. The assembly of nodal KCL equations is written in matrix form as
3.6where ** I** is the vector of external input currents,

**is the vector of nodal voltages and**

*V***is the global coefficient matrix. The global coefficient matrix is built from each resistor element in the network using connectivity information. After assembling the global coefficient matrix, the essential boundary conditions are applied in order to make the system unique and non-singular. For this model, the essential boundary conditions are user-defined voltage values and are applied to the top and bottom of the RVE in the traditional FEA manner. The system is then solved for the unknown electric potentials at the CNT nodes. These voltages,**

*K**V*

_{i}and

*V*

_{j}, are then post-processed to determine the current in each linear resistor,

*R*

^{e}, using equation (3.3). This solution technique is common in FEA and can be found in standard FEM books. Element currents provide the basis for two important results presented by the model. The first result relates to backbone identification. Not all resistor elements, or fibres, in the network contribute to electrical properties. In fact, only a small fraction will make up the path of least resistance and carry the bulk of the electrical current. This fraction of nanotubes is known as the backbone [45]. Theoretically, the backbone consists of all fibres with non-zero currents [53]. Besides identifying the backbone, the model simultaneously calculates the effective electrical properties of the composite. The total current entering the RVE is calculated by summing the element currents of resistors in contact with the top of the microstructure. Finally, the effective electrical properties of the system are calculated using the equations [51,52] 3.7 3.8 3.9

Equation (3.7) yields the effective electrical conductance of the system in units of siemens (S), where *V* _{top} and *V* _{bottom} are the applied voltages on the RVE boundaries. Conductance measures how easily current flows in the system, and involves the inverse of Ohm's law in equation (3.8). Resistance measures strength of opposition to current flow in the system and is analogous to friction in mechanical systems. Equation (3.9) gives the effective electrical conductivity of the nanocomposite, where *t* is the RVE thickness.

Clearly, the approach has a large number of user-specified input parameters; this enhances its versatility and flexibility. Again, as in the preceding case of mechanical and thermal properties, a Monte Carlo approach with adequate sampling is used for a reliable determination of the electrical properties.

### (e) Numerical results

The numerical model contains six distinct sources of stochastic behaviour. Each is a random variable either sampled from a distribution or directly calculated from another stochastic variable. These inputs are readily controllable for the particular application or study. A study on the stochastic behaviour of this model determined the impact of the random CNT length and diameter variables. First, CNT length and diameter are sampled from Weibull and lognormal distributions, respectively. For the second set, CNT dimensions are fixed at 149.5 and 1.77 nm for length and diameter, respectively. These values maintain the same average aspect ratio for both sets of data. This allows the statistical variation between the two sets to be restricted to the variables themselves and to isolate the effects they have on the results. The impact of fixing CNT dimensions on average effective conductivity is shown in figure 13.

Results in figure 13 confirm that eliminating the randomness in CNT dimensions has minimal impact on the percolation ratio and electrical conductivity of the nanocomposite. Both sets percolate all of the 500 RVEs in the volume fraction range of 0.008–0.024, with most data points differing by a few per cent from each other. The effective conductivity is nearly the same, with only one data point that differs by an order of magnitude. The incorporation of all six stochastic features of the model would theoretically produce the most accurate representation of the real-world nanocomposite system. For the purpose of this study, however, it is only important to demonstrate the derivation of reliable qualitative results and of quantitative estimates of effective electrical properties. Thus, it is beneficial to reduce the number of random variables in the model to reduce the computational cost of the program.

Next, a study is conducted on the phenomenological nature of atomic-scale effects on electron transport. The model incorporates this physical attribute with the input variable of tunnelling distance. Four sets of simulations are conducted with different tunnelling distances ranging from 0 to 10 nm. Comparisons of the percolation ratio and electrical conductivity results are presented in figure 14.

It is readily seen from figure 14 that the percolation ratio and effective conductivity results mirror the physical behaviour of the material. As the tunnelling distance increases, the ratio of percolated RVEs increases for each volume fraction. Increasing the tunnelling distance allows more CNTs that previously had not formed connections to meet the bonding criterion. Thus, the probability of spanning network formation at a given volume fraction increases. The effective conductivity also increases at each volume fraction, as shown in figure 14*b*. This result reflects an increase in RVEs containing complete spanning networks and an increase in the number of percolation paths available for current to flow. The differences in percolation and conductivity results are amplified for the case of 10 nm. The significantly larger tunnelling distance causes a pronounced reduction in the percolation threshold. The material is completely insulating at 0.008 volume fraction when tunnelling is less than 10 nm. However, at 10 nm, 1 per cent of the RVEs form complete networks at the same concentration. Lower percolation thresholds and higher effective electrical conductivities are fundamental aspects of the electron tunnelling phenomenon in nanocomposites.

The electric potential and flux can be visualized in the RVE. For this study, CNTs are assumed to be wavy with fixed lengths and diameters. They are randomly dispersed in a 1 μm^{2} RVE at a volume fraction that is high enough for complete percolating networks to form. The first step to determine nodal voltages and element fluxes is to identify the complete spanning network of CNTs. The spanning network algorithm identifies all CNTs in contact or within tunnelling distance of each other, and establishes connection points for the creation of a resistor network. The percolation network identified for two random RVEs is presented in figure 15.

Figure 15 shows that there are continuous CNT structures spanning the longitudinal direction of both RVEs. Each highlighted CNT meets the bonding criteria with either another nearby CNT or the voltage source at the top of the RVE. Consequently, there are fibres included with the spanning network that clearly will not carry any current and may not be near another fibre. With the spanning networks in hand, the equivalent resistor network is solved. The nodal voltages are presented in figure 16.

The results in figure 16 show that RVE no. 1 has more CNTs at a high voltage than RVE no. 2. Many have potentials close to 100 V that extend more than half way to the bottom of the RVE. Also, the transition from source to drain voltage is much smoother for RVE no. 1, with several CNTs in the mid-voltage range. RVE no. 2 is characterized by a very sharp voltage transition near the source. These transitions, or cliffs, appear to occur in regions with few connecting CNTs [55]. A key feature of the model is its ability to identify the conducting backbone of a complete CNT network. Simply stated, the backbone is the current-carrying part of the percolating network [45]. The final step in this process is to identify this backbone. Results show that the largest currents in RVE no. 1 occur close to the left side, and in RVE no. 2 they are towards the right side, as depicted in figure 17.

The backbone can be readily identified in each RVE in figure 17. The regions of relatively high flux indicate the path of least resistance through the entangled mass of CNT fibres. Current drops are visible in the backbone. These indicate the existence of parallelism in the circuit where charge flows through multiple CNTs at the same time near a junction. Though each RVE in this study contains one backbone, it is possible for multiple parallel backbones or branches to occur at higher volume fractions.

The percolation probability results shown are for a typical multi-walled carbon nanotube/ polycarbonate (MWCNT/PC) material. Inputs to the model reflect those of real constituent materials and are summarized in table 1.

The RVEs are sampled at volume fractions in the range 0.001–0.0045 and the CNTs are wavy. Besides reporting the percolation probability at several volume fractions, a method for accurately interpolating the results at all volume fractions is introduced. In this context, note that Li & Chou [86] proposed that the percolation probability of a composite with random arbitrary filler shapes can be approximated by the cumulative distribution function (CDF) of the volume fraction according to the equation
3.10In this equation, *f* is the volume fraction, *μ* is the mean percolation volume fraction and *σ* is the standard deviation (s.d.). The method used by Li & Chou [86] to approximate probability using a CDF is implemented in this model to provide information at all volume fractions. The mean and s.d. of the percolation volume fraction are 0.001 8128 and 0.000 195 327, respectively. Using equation (3.10), the results are fitted with the CDF for a normal distribution and are shown in figure 18.

The CDF fit shown in figure 18 is in good agreement with the model results in the percolation region. The region spans volume fractions of 0.0012–0.0026, with both the data and CDF curves converging to 1.0 by 0.0026 volume fraction. The volume fractions at which probabilities occur differ by up to 0.0001 between the model and the CDF fit for this material. However, at this low CNT concentration, the difference is a matter of a few MWCNTs. Although there are clear differences between the data and CDF fit, the method is a simple and useful estimate of the probability of percolating at a given CNT concentration for a simulated material.

The effective electrical conductivity results include comparison of data with experiments and presentation of a new approach to characterize CNT-based nanocomposites. The electrical conductivity results appear to follow a power-law curve. According to the percolation theory, the conductivity in the threshold region can be fitted by the percolation power-law equation [83,87]
3.11The variable *p* is CNT volume fraction, *p*_{c} is the percolation threshold, *σ*_{o} is a coefficient that depends on nanotube conductivity, and *t* is the critical exponent. Equation (3.11) expresses the dependence of a macroscopic property, effective conductivity, on microscopic properties, namely volume fraction and percolation threshold. The relationship between them is characterized by the critical exponent, *t*. The critical exponent is determined from the reduced volume fraction plot shown in figure 19.

The critical exponent and coefficient are 1.702 and 4.7×10^{4}, respectively. Figure 19*a* shows a good power-law fit with an *R*^{2} value of 0.9964. The resulting power-law fits very well with the original simulation data points for volume fractions above the percolation threshold, *p*_{c}.

Effective conductivity data from several experiments are compared with the numerical data. Experiments chosen for effective conductivity comparison include Ramasubramaniam *et al.* [88], Ounaies *et al.* [82], Skákalová *et al.* [16], Hu *et al.* [89], Sandler *et al.* [80] and Gojny *et al.* [90]. These experiments provided a wide range of nanocomposite combinations for comparison. All types of CNTs were included in combination with five different matrix materials. Further, several different processing methods were represented by the data. The coordinate finding program [42] is used to gather data points for the experimental comparisons. The model proposed herein simulates 500 random RVEs of wavy MWCNTs embedded in a PC matrix material. A 100 V potential difference is applied across the boundaries and the contact resistances are sampled stochastically. Comparison of the effective conductivity results is shown in figure 20.

Figure 20 illustrates the wide variation in the reported experiments. The model results appear to be close in magnitude to the experiments. Unfortunately, few data points exist for each material. This makes a close comparison difficult. Regardless, it is important to note that the model follows the same trend as the experimental results. Both have the characteristic jump in effective conductivity near the percolation threshold and continue to increase modestly according to the percolation power law. Because CNT type and morphology dominate the effective electrical properties, it is important to closely examine material results using MWCNTs. Only the data of Hu *et al.* [89] in figure 20*a* give a percolation threshold close to the model. The MWCNT materials of Sandler *et al.* [80] and Gojny *et al.* [90] in figure 20 all have thresholds at lower volume fractions. Percolation threshold and critical exponent are compared with several reported experiments in figure 21.

Figure 21 shows that results from the model are in good agreement with experiments. Experimental percolation thresholds are between 0.008 and 0.028 volume fraction. The model estimates a percolation threshold of 0.001 17 for a CNT aspect ratio of 500. The relatively high-aspect-ratio nanotubes used in these composites percolate at very low volume fractions. The critical exponent, *t*, for the same simulation set overestimates the critical exponent. It is predicted to be 3.63, whereas the experimental values range from 1.3 to 3.3. Although the predicted critical exponent seems high, it is within the 1.3–4 range frequently reported in experiments [104]. Furthermore, the values tend to be dimensionally dependent, with *t*=1.33 for two-dimensional and *t*=2 for three-dimensional models [104]. Therefore, the high value from the pseudo-three-dimensional model suggests its capability to reflect the electrical behaviour of a full three-dimensional model.

Benefits of the current model include the capacity to estimate percolation probability and the effective electrical conductivity at volume fractions in the percolation region. The CDF fitted to percolation ratio data predicts the percolation probability using equation (3.10). Therefore, the likelihood that the input material will form conductive networks of CNTs at any given volume fraction is known. Further, the effective electrical conductivity is estimated from the power-law fit according to equation (3.11). By combining these results-based fits, the electrical characteristics of the input material are known at any volume fraction in the percolation region. The first material characterized is a typical SWCNT/epoxy or SWCNT/alumina with filler aspect ratio of 100. Table 2 summarizes this input nanocomposite.

The model predicts a percolation threshold of 0.009 volume fraction, a critical exponent of 2.4834 and a coefficient of 6.59×10^{6}. The CDF fit has a calculated mean volume fraction of 0.017 and s.d. of 0.003 11. The fits are interpolated over the percolation threshold region of 0.008–0.03 volume fraction. Materials for electrostatic discharge (ESD) applications typically require conductivities of at least 10^{−6} S m^{−1} [102,106,107]. At these low conductivities, the material gradually dissipates electrical charge build-up on devices such as electronic components. For applications involving electromagnetic interference (EMI) shielding, conductivities greater than 1 S m^{−1} are necessary [88]. Electrical conductivity and percolation probability are visualized in figure 22.

Figure 22 shows a measure of a nanomaterial's electrical performance. At low volume fractions, such as 0.01, the given SWCNT/epoxy material is expected to have an effective conductivity of about 0.24 S m^{−1}. Therefore, this material is suitable for ESD applications. However, it has only a 1.22 per cent probability of forming the spanning CNT networks required to achieve that conductivity. The probability that conducting networks of CNTs exist increases rapidly over a small range of volume fractions. Comparatively, the effective electrical conductivity of the materials with spanning networks increases from the ESD to EMI range over the same range of volume fractions.

Finally, a MWCNT/PC material is simulated with 500 random RVEs. The material's electrical properties are predicted over a volume fraction range of 0.001–0.0045. This material has a much larger aspect ratio and lower CNT conductivity compared to the SWCNT/epoxy material analysed previously. Percolation probability has a CDF fit with calculated mean volume fraction of 0.001 81 and s.d. of 0.000 195. A percolation threshold of 0.001 17 volume fraction, a critical exponent of 3.63 and a coefficient of 4.6×10^{8} are calculated for the power-law fitted curve. Both fits are interpolated over the volume fraction range in figure 23.

The electrical conductivity for this material is predicted to never exceed the required amount for EMI applications at volume fractions below 0.0045. Rather, it is obvious from figure 23 that this material is best suited for ESD applications. Also, the material percolates at a volume fraction nearly 10 times lower than the SWCNT/epoxy nanocomposite. The lower effective conductivity and percolation threshold are probably due to the lower intrinsic conductivity and higher aspect ratio of MWCNTs over SWCNTs.

## 4. Concluding remarks

This paper has discussed a model to predict the nonlinear elastic and thermal properties of SWCNT-reinforced composites using the EFFEM. An accurate stress–strain curve for these nanocomposites can be generated efficiently by incorporating the nonlinear modulus of elasticity of the polymer matrix using an incremental approach. The accuracy of the numerical results has been assessed by comparison with experimental results available in the literature. The results have been found to be accurate within 5 per cent for small strains less than 5 per cent. The effective conductivity of the composite has also been investigated while incorporating the nonlinear thermal conductivity of SWCNTs. Using an iterative approach, the numerical results have been compared with several experiments in the literature. The approach has been shown to be quite accurate for small volume fractions—less than 3 per cent—but inaccurate for large volume fractions. This has been attributed to complex mechanisms inhibiting heat conduction at the nanotube–matrix interface.

The paper has also addressed the determination of the effective electrical properties of SWCNT-reinforced composites. For this purpose, a spanning network algorithm and the resistor network method have been used. The spanning network algorithm has incorporated the electron tunnelling effect to identify complete CNT networks. The results of this algorithm have led to the estimation of the percolation probability in the material. Finally, the resistor network method has been used to represent the CNTs and contact points in the identified network. KCL and the FEM have been used to assemble and solve the system equations. From the results, the current-carrying backbone and effective conductivity of the RVE have been calculated at each volume fraction. The electrical model has used the same method as the elastic and thermal model to develop the RVE morphology, but has differed in how the problem was discretized. When determining elastic and thermal properties, the model used FEM to add the physical properties to the polymer matrix. Conversely, when determining the electrical properties, the model used an equivalent resistor network and FEA to solve for the unknowns. It is important to note the need for different discretization processes for determining different physical properties. The model has produced effective electrical conductivities and percolation probabilities for a variety of CNT-based nanocomposites in the percolation region. The predicted effective electrical conductivities, percolation thresholds and critical exponents for several specific materials have compared well with several experiments. Therefore, the accuracy and versatility of the model have been shown. The accuracy of the model, however, could be improved with better understanding of both the nanoscale physics that governs electron transport and the impact that manufacturing techniques have on electrical properties.

Monte Carlo simulations have minimized the statistical variation in the derived numerical results both for the determination of the effective mechanical and thermal properties and for the determination of the effective electrical properties of the nanocomposite.

## Acknowledgements

The partial financial support of this work through a contract from the Department of Defense to Clarkson Aerospace–Minority Leaders Programme is acknowledged with pleasure.

## Footnotes

One contribution of 17 to a Theme Issue ‘A celebration of mechanics: from nano to macro’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.