We report Monte Carlo simulations of systems of polydisperse prolate and oblate ellipsoids using the critical path based tunneling-percolation model. For polydisperse prolate ellipsoids, the critical percolation volume fraction ϕc is shown to have a quasi-universal dependence on weight-averaged aspect ratio. For polydisperse oblate ellipsoids, ϕc is shown to have a quasi-universal dependence on the apparent aspect ratio, which is a function of up to fourth moment of the size distribution, as given by the generalized connectedness percolation theory. The functions are observed to approach the theoretical predictions for higher volume fractions and higher aspect ratios. The model predictions are compared with experimental data available on polydisperse multi-walled nanotubes (prolate ellipsoids) and graphene nanoplatelets (oblate ellipsoids) to estimate the tunneling lengthscale which is found to be well within the expected range.

Nanofillers with highly anisotropic shapes, such as carbon nanotubes (CNTs), graphene nanoplatelets (GNPs), carbon black (CB), and metallic nanowires,1 are used as inclusions in polymer matrix materials to generate nanocomposites with superior electrical, mechanical, and thermal properties. In this study, the focus is on the enhanced effective electrical conductivity. Electrically conductive polymer nanocomposites can be used in static dissipative, semi-conductive, conductive, or electromagnetic shielding applications. The increased conductivity is also desired in preventing the aircraft structures made of carbon based polymer composites from localized damage due to lightning strikes.2 The strong dependence of the electrical conductivity on interfiller separation has also been exploited by using nanocomposites for strain measurement applications.3–5 Given the wide range of applications, electrical conductivity of nanocomposites has been an intense topic of research over the past few decades as an in-depth understanding of the various factors governing the electrical behavior is a crucial step in the advancement of nanocomposites based technology.

Electrical conductivity of insulating polymers can be increased by several orders of magnitude by adding a small amount of highly conductive nanofillers. The sharp rise in the electrical conductivity of polymer nanocomposites results from formation of a percolating network of electrically interconnected nanofillers through electron tunneling.6,7 For a given volume fraction of nanofillers, the polymer matrix layer separating two fillers is sufficiently thin so as to allow electrons to tunnel across the barrier. The tunneling conductance between two nanoparticles can be expressed as8 

gij=g0exp(2δijξ),
(1)

where δij is the shortest distance between two particles i and j, ξ is the tunneling lengthscale, and g0 is a constant proportional to the contact resistance.

Given the extreme contrast between the electrical conductivity of the fillers and the polymer matrix, the effective electrical conductivity of the nanocomposite can be modeled as a tunneling-percolation problem with conductive fillers distributed in a continuous non-conductive matrix. Modeling this behavior is a challenging task given the large number of variables (i.e., “all thinkable parameters” as concluded in Ref. 9) that play a key role in determining the percolation behavior of the system.

It was shown in Ref. 8 that effective electrical conductivity of nanocomposites can be obtained from the global tunneling network (GTN) of the interparticle tunneling conductors as defined in Eq. (1). Furthermore, as the tunneling conductance varies exponentially with the interparticle distance, the effective conductivity of the GTN can also be obtained using the critical path approximation.10 The major advantage of the critical path method is that the complexity of the problem of finding the effective conductivity of the network is significantly reduced. Instead of finding the effective conductivity by numerically decimating the GTN, one can estimate the conductivity geometrically using

σσ0exp(2δcξ),
(2)

where δ/2 is the thickness of the soft penetrable shell surrounding the hard particles and δc is the smallest δ for which the fillers at a given volume fraction form a percolating network. A model based on the critical path method,8 as explained in detail in Section II, is used in this study.

Aspect ratio of nanofillers is identified as the primary factor that has a significant impact on the electrical conductivity of conductive polymer nanocomposites.1 The effect of aspect ratio on the conductivity has been extensively studied through experiments,11 modeling,8 and theoretical investigations.12,13 However, in these studies, it is often assumed that all the fillers are identical in size and shape to simplify the analysis.

Comparison of the simulations or theoretical predictions of the idealized models with experimental data is not straightforward. Many other factors such as the size distribution of fillers, changes in shape, size, and aspect ratio due to breaking of fillers during processing, and agglomeration of fillers lead to a system that differs significantly from the idealized model8 with uniformly dispersed identical spheroidal fillers.1 As a result, the idealized model may fail to represent the actual composite and deviation from experimental data can be observed. Therefore, understanding the effect of such parameters is highly important when comparing the results of the model with experiments, and more so, when using the model as a predictive tool.

Analytical studies on polydisperse rods, based on the mean-field assumption,9,13–15 have shown that the percolation threshold ϕc is extremely sensitive to polydispersity. It was found that ϕc is inversely proportional to the weight average of the length distribution Lw=L2/L, for penetrable and impenetrable systems of rods. Therefore, with increasing polydispersity, the ϕc decreases even if the mean length L is fixed. The longer rods contribute strongly in establishing the percolating network, especially when the length difference between longer and shorter rods increases. The analytical models, although valid strictly in the mean-field sense and in the slender-rod limit (high aspect ratios), provide a powerful way to interpret the experimental and simulation results in a polydispersed system of elongated fillers like CNTs.

The effect of length polydispersity on a system of spherocylinders with a fixed diameter D is studied in Ref. 16 using Monte Carlo (MC) simulations. Spherocylinders are cylinders with diameter D, length L, and hemispherical caps of diameter D at the ends. For a polydispersed system of penetrable (soft-core) spherocylinders, ϕc is shown to be inversely proportional to Lw, as predicted theoretically in the slender-rod limit.9,13–15 For a sufficiently large L2/D, a universal dependence of the form ρδc3=F(L2/δc) has been observed, where ρ is the number density. Mutiso et al.17 have used a generalized form of the excluded volume model and percolation simulations of soft-core polydispersed cylinders to model the effect of finite L/D and size polydispersity on ϕc. Their model is not limited by the slender-rod limit14 and works for finite, but penetrable, rods as well. Chen et al.18 have presented a generalization of the excluded volume model12 that can predict the percolation threshold for penetrable polydisperse spherocylinders. Recently, Meyer et al.19 studied the effect of polydispersity on percolation behavior of rod-like inclusions using theoretical and Monte Carlo simulation approaches. They modeled Inclusions as spherocylinders with a hard core and a penetrable shell and accounted for polydispersity in shell thickness of the fillers as well as length and diameter. Such analytical and simulation based studies have greatly enhanced our understanding of percolation behavior of penetrable and impenetrable polydisperse rodlike fillers.

Based on a generalized connectedness percolation theory, Otten and van der Shoot9 have theoretically predicted the percolation threshold for hard polydisperse platelike objects with soft penetrable shells. For monodisperse hard discs, ϕc is shown to be independent of the aspect ratio. However, their model predicts that ϕc of polydisperse hard discs depends on up to fourth moment of the size distribution, and so a strong influence of the diameter polydispersity is expected. However, the theoretical prediction is only valid in the limiting case of thin discs.

In this study, we analyze the effect of size polydispersity on percolation of randomly oriented prolate and oblate ellipsoids for a wide range of aspect ratios using MC simulations and compare our results with the theoretical predictions. The results for polydisperse prolate ellipsoids are shown to be in agreement with the theoretical predictions9 and polydisperse spherocylinders simulations.16 We show that the results for monodisperse and polydisperse oblate ellipsoids can be collapsed on a universal functional form when the aspect ratio is sufficiently large, but still much smaller than the thin plate limit. Therefore, the theoretical predictions based on connectedness percolation9 can be confidently used in practice where the aspect ratios of the fillers used is sufficiently large.

The paper is organized as follows. The tunneling-percolation model is presented in Section II. The MC results for polydisperse prolate and oblate ellipsoids representing CNTs and GNPs, respectively, are compared with the theoretical model in Section III, and with experimental results in Section IV, followed by conclusion in Section V.

The fillers are modeled as impenetrable spheroids with aspect ratio α that best captures the morphology of the filler. α is defined as the ratio of polar to equitorial semi-axes of the spheroid. For CNTs, the cylinder-like shape is obtained using α1, while for GNPs the plate-like shape is obtained with α1. The random distribution of fillers is generated using the random sequential addition (RSA) method.20 In this method, the desired fillers are added one by one such that the recently added particle does not collide with the previously added ones. The test of collision between two ellipsoids is implemented using a procedure proposed by Perram and Wertheim.21 The process is continued until the desired volume fraction is reached. Although simple, RSA suffers from a major drawback that for a given α there is a theoretical limit on the maximum volume fraction that can be attained by finite random trials.20 However, the maximum attainable volume fractions by RSA were found to be sufficient for the current study. It is also worth noting that this technique generates configurations that are not equilibrated, which can be problematic for configurations with high density of fillers. Nevertheless, the discrepancy is not significant at small volume fractions22 that are considered in the present work. MC relaxation technique is generally used to obtain equilibrated configurations. Also, percolation behavior of RSA generated monodisperse oblate ellipsoids with and without MC relaxations was found not to differ significantly in Ref. 8.

The position of the filler is chosen within a cube using three random numbers between 0 and l, where l is the edge length of the cube. Periodic boundary conditions are enforced along all directions. The random orientation of the polar axis of spheroids is generated with a uniformly distributed random variable on a unit sphere surface using θ=2πm and ϕ=arccos(2n1), where m and n are random numbers in [0,1], while θ and ϕ are the spherical coordinates giving the polar and azimuthal angles of a point on the unit-sphere surface, respectively.

In the simulations performed, the box size l is chosen such that it is at least 5 times the largest filler dimension and that the total number of fillers is at least 25 000. In this way, it is ensured that the finite-size effects does not significantly affect the percolation based results. The box is further divided into cubical sub-blocks of the size of the largest filler dimension. Such a subdivision narrows down the collision check only to the surrounding sub-blocks. The expensive ellipsoid collision check is bypassed if the distance between the centers of the two fillers is larger than the sum of their largest semi-axes, thus further improving the speed and efficiency of the collision check operation.

After generating the RSA configuration of fillers with a desired volume fraction ϕ, all fillers are coated with a soft penetrable shell of uniform thickness δ/2. The objective is to find δc which is the minimum δ for which a percolating cluster of the interconnected fillers spanning the entire domain is formed. Here, for simplicity, the underlying assumption is that the tunneling mechanism is isotropic, i.e., has no directional dependence.

The percolation check is performed using the tree-based union/find algorithm with path compression proposed by Newmann and Ziff.23 In this method, the hard core fillers with a penetrable shell of thickness δ/2 are added one by one to the box. After adding each filler, the cluster information, stored using a tree datastructure, is updated based on the contacts established with the neighboring fillers through the penetrable soft shell. The objective is to check if for a given ϕc and δ, a percolating cluster spanning the cube in a given direction is formed. The details of the percolation check algorithm and its implementation can be found in Ref. 23.

Starting with an initial guess for δ, for a RSA generated filler configuration at a desired volume fraction ϕ, the percolation check is performed as described above. If the system is found to percolate, then δ is decreased, and if percolation is not achieved, then δ is increased. Therefore, δc is obtained to within a desired accuracy using a bounded root finding methods such as the bisection method.

The nanofiller mixing processes such as sonication and shear mixing, which improve the filler dispersion, lead to significant changes in the size distribution of the CNTs Refs. 7 and 24. Experimental studies that have analyzed CNT size distribution have shown that Weibull distribution effectively captures the polydispersity.25,26 Such distributions with significantly higher moments have a notable effect on the percolation behavior, and hence the effective conductivity of the composite.

As mentioned before, the effect of polydispersity on rodlike fillers has been investigated using analytical methods9 and MC simulations.1,16,17 In this section, the results for polydisperse prolate ellipsoids are compared with the theoretical predictions in the slender-rod limit. The analysis is along the same lines as that of MC simulations of spherocylinders in Ref. 16, and we expect to see similar trends for prolate ellipsoids, especially at higher aspect ratios.

Using the liquid-state integral equation theory, the critical percolation volume fraction ϕc for polydisperse impenetrable rods with diameter D and length L is given as9 

ϕcπ4ρD2L=D2L2δcL2.
(3)

Then,

ρδc3=2δc2πL2.
(4)

The theoretical prediction, applicable in the slender-rod limit, i.e., L/D1, was used in Ref. 16 to analyze the percolation behavior of impenetrable spherocylinders. It was shown that, for different polydisperse systems with the same L2/D, the ρδ3 values are essentially independent of the actual size distribution. Moreover, with increasing L2/D, the simulation data tend to converge towards a common functional form which approaches Eq. (4) for L2/δc1.

Similar arguments can be extended for monodisperse and polydisperse systems of prolate fillers studied here (Fig. 1). The semi-axes of the prolate ellipsoids are given as {L/2,D/2,D/2}. Prolate ellipsoids following a bimodal distribution f(L)=pδ(LL1)+(1p)δ(LL2) are used to account for length polydispersity, where p is the number fraction. All prolate ellipsoids are assumed to have the same D, and all lengths are normalized with D. The aspect ratio α is given as L / D. Similar to spherocylinders,16 different size distributions with the same L2/D collapse on a single curve as shown in Fig. 2. Moreover, for larger L2/D, the collapse of data suggests a common functional form that is independent of L2/D. Therefore, as expected, the quasi-universal behavior for a given L2/D that was reported for spherocylinders16 is observed for polydisperse prolate ellipsoids as well. A similar analysis on polydisperse oblate ellipsoids is presented next.

FIG. 1.

Polydisperse prolate ellipsoids following a bimodal distribution with α1=10,α2=20, and p = 0.5 at ϕ=0.025. The fillers shown in green belong to the percolating cluster.

FIG. 1.

Polydisperse prolate ellipsoids following a bimodal distribution with α1=10,α2=20, and p = 0.5 at ϕ=0.025. The fillers shown in green belong to the percolating cluster.

Close modal
FIG. 2.

ρδc3 as a function of L2/δc for monodisperse and polydisperse prolate ellipsoids. The data having the same L2/D can be observed to collapse well. The theoretical prediction for the slender-rod limiting case Eq. (4) is shown as a dotted line for comparison.

FIG. 2.

ρδc3 as a function of L2/δc for monodisperse and polydisperse prolate ellipsoids. The data having the same L2/D can be observed to collapse well. The theoretical prediction for the slender-rod limiting case Eq. (4) is shown as a dotted line for comparison.

Close modal

The inverse proportionality between ϕc and aspect ratio that is observed for prolate ellipsoids is absent in the case of oblate ellipsoids.9,27,28 After an initial drop in ϕc with increasing aspect ratio, ϕc saturates to a limiting value independent of the aspect ratios. This behavior of a system of disc-like fillers is suggested to be the outcome of interference of the percolation transition with the isotropic-nematic phase transition.28 However, ϕc is shown to be significantly altered in the presence of diameter and thickness polydispersity.9 

The analytical expression for ϕc for polydisperse hard discs based on connectedness percolation theory and second-virial approximation to the pair connectedness function was derived in Ref. 9 as

ϕc=4LD2δc[ 1(5π+6)D2+16π(π+6)DD3+(π+6)2D4+8π(π6)D22 ]=4LD2δcA,
(5)

where L is the disc thickness and D is the disc diameter. Higher moments of the size-distribution appear in the expression, unlike in the slender-rod case, where only the first and second moments appear. Equation (5) is applicable when δcD and when L and δc are of the same order of magnitude. Unlike for rodlike particles, the second-virial approximation is supposed to break down for platelike and spherical particles as the contribution of the higher order terms cannot be neglected.9 However, for δD, it turns out to be a reasonable approximation for spheres, and hence is presumed to be applicable to platelike particles as well when deriving Eq. (5).9 Therefore, Eq. (5) is expected to work at least qualitatively within the second-virial approximation. Moreover, ϕc predicted using Eq. (5) was shown to match well with experimental data29 where the GNP size polydispersity information was obtained using dynamic light scattering.

For a polydisperse distribution of hard discs, ϕc=ρπD2L/4. Then, using Eq. (5)

ρδc3=16πδc2A,
(6)

where 1/A represents the expression in the square brackets of Eq. (5). Therefore, in the thin disk limit, scaling of the form Eq. (6) is expected.

Using the percolation-tunneling model of ellipsoidal fillers, we have simulated monodispersed and polydispersed configurations of oblate ellipsoids (Fig. 3). The semi-axes of oblate ellipsoids are represented as {D/2,D/2,L/2}, and the aspect ratio is given as α=L/D. The results are presented by plotting ρδc3 as a function of A/δc. The plots for monodisperse cases with different aspect ratios are shown in Fig. 4. With increasing aspect ratio, the curves are observed to converge on to a universal functional form ρδc3=F(A/δc). Equation (6) for the limiting case of thin hard discs is also presented for comparison. The simulated oblate ellipsoids approach the thin disc behavior when α is very small and δc/D1.

FIG. 3.

Polydisperse oblate ellipsoids following a bimodal distribution with α1=1/10,α2=1/20, and p = 0.5 at ϕ=0.025. The fillers shown in green belong to the percolating cluster.

FIG. 3.

Polydisperse oblate ellipsoids following a bimodal distribution with α1=1/10,α2=1/20, and p = 0.5 at ϕ=0.025. The fillers shown in green belong to the percolating cluster.

Close modal
FIG. 4.

ρδc3 as a function of A/δc for monodisperse oblate ellipsoids with aspect ratios α={1/2,1/10,1/20,1/40,1/100}. For the monodisperse case, A=(10π+12)D. The analytical expression for the thin hard discs, Eq. (6) is also shown as a dashed line for comparison.

FIG. 4.

ρδc3 as a function of A/δc for monodisperse oblate ellipsoids with aspect ratios α={1/2,1/10,1/20,1/40,1/100}. For the monodisperse case, A=(10π+12)D. The analytical expression for the thin hard discs, Eq. (6) is also shown as a dashed line for comparison.

Close modal

A bidisperse system of oblate ellipsoids described by a bimodal distribution f(D)=pδ(DD1)+(1p)δ(DD2) is simulated, where D1>D2 and p[0,1] is the number fraction of the larger ellipsoids. To simplify the analysis and to focus on polydispersity in the larger dimension, all oblate ellipsoids in the simulations have the same L, i.e., all GNPs are assumed to have the same thickness. Also, all dimensions are normalized with L.

The effect of polydispersity on the critical distance δc can be examined by using the reduction factor R, which is the ratio of δc to δc0, where δc0 is the critical distance for the monodisperse system of ellipsoidal fillers with semi-axes equal to the mean semi-axes of the polydisperse case, i.e., {D/2,D/2,L/2}. For the bimodal distribution, D=pD1+(1p)D2.

An analytical expression for the reduction factor R0 is obtained using Eq. (5) as

R0=δcδc0=2(5π+6)D2A,
(7)

where

δc0=2Lϕc(5π+6)
(8)

is the critical distance for monodisperse hard discs at a volume fraction ϕc. Both δc and δc0 are for the same volume fraction ϕc.

The reduction factors R obtained from simulating oblate ellipsoids following the bimodal distribution for D2/L=5 and D1/D2=4 are plotted in Fig. 5 as a function of p. The plots are generated for different total volume fraction ϕc of the fillers. The maximum ϕc used is limited by the RSA procedure used to generate the realizations. It can be observed that R obtained from the simulations becomes qualitatively similar to R0 as ϕc increases. Similar trends can be observed in Fig. 6 for D2/L=5 and D1/D2=2. While for small ϕc, the condition necessary for the theoretical prediction R0 to hold, i.e., δcD is not satisfied, and hence there is a large difference between R and R0. With increasing ϕc,δc/D decreases and the theoretical result is approached by the simulated data. Apart from ϕc, the aspect ratio D2/L is also expected to affect R. The effect of D2/L on R is shown in Fig. 7, where R is plotted as a function of p for ϕc=0.025,D1/D2=2, and D2/L={2,4,10,20}. It is again observed that with increasing D2/L, R approaches the analytical prediction R0.

FIG. 5.

The reduction factor R plotted as a function of the number fraction p of larger oblate ellipsoids from a bimodal distribution with D2/L=5 and D1/D2=4 for increasing ϕc. The theoretically predicted reduction fraction R0 according to Eq. (7) is shown as a solid line.

FIG. 5.

The reduction factor R plotted as a function of the number fraction p of larger oblate ellipsoids from a bimodal distribution with D2/L=5 and D1/D2=4 for increasing ϕc. The theoretically predicted reduction fraction R0 according to Eq. (7) is shown as a solid line.

Close modal
FIG. 6.

The reduction factor R plotted as a function of the number fraction p of larger oblate ellipsoids from a bimodal distribution with D2/L=5 and D1/D2=2 for increasing ϕc. The theoretically predicted reduction fraction R0 according to Eq. (7) is shown as a solid line.

FIG. 6.

The reduction factor R plotted as a function of the number fraction p of larger oblate ellipsoids from a bimodal distribution with D2/L=5 and D1/D2=2 for increasing ϕc. The theoretically predicted reduction fraction R0 according to Eq. (7) is shown as a solid line.

Close modal
FIG. 7.

The reduction factor R plotted as a function of the number fraction p of larger oblate ellipsoids from a bimodal distribution with D1/D2=2 at ϕc=0.025 and D2/L={2,4,10,20}. The theoretically predicted reduction fraction R0 according to Eq. (7) is shown as a solid line.

FIG. 7.

The reduction factor R plotted as a function of the number fraction p of larger oblate ellipsoids from a bimodal distribution with D1/D2=2 at ϕc=0.025 and D2/L={2,4,10,20}. The theoretically predicted reduction fraction R0 according to Eq. (7) is shown as a solid line.

Close modal

The data for oblate ellipsoids following a bimodal distribution with α1=L/D1=1/10,α2=L/D2=1/2, and p={0,0.5,1} obtained for a range of ϕc is shown in Fig. 8, where ρδc3 is plotted as a function of A/δc. The data for monodisperse oblate ellipsoids with D/L=A/(10π+12)/L, where A is calculated for the bimodal distribution given above at p = 0.5 is also shown. For α1=L/D1=1/10,α2=L/D2=1/2, and p = 0.5, D/L=7.435. For a given A/L, the ρδc3 against A/δc curves collapse well and can be concluded to be independent of the specific size-distribution. The data for a wide range of monodisperse and bidisperse oblate ellipsoids simulated in this study are collectively shown in Fig. 9. For oblate ellipsoids with sufficiently large A/L, the curves can be observed to collapse well, suggesting a universal functional form ρδc3=F(A/δc). Also, as Aδc and for a system of hard discs, the scaling function approaches Eq. (6).

FIG. 8.

ρδc3 as a function of A/δc for oblate ellipsoids following a bimodal distribution with α1=L/D1=1/10,α2=L/D2=1/2, and p={0,0.5,1}. The data for p = 0.5 collapses well with the data for a monodisperse case having the same A/L value.

FIG. 8.

ρδc3 as a function of A/δc for oblate ellipsoids following a bimodal distribution with α1=L/D1=1/10,α2=L/D2=1/2, and p={0,0.5,1}. The data for p = 0.5 collapses well with the data for a monodisperse case having the same A/L value.

Close modal
FIG. 9.

ρδc3 as a function of A/δc for monodisperse and bidisperse oblate ellipsoids. The theoretical prediction according to Eq. (6) is shown as a dotted line for comparison.

FIG. 9.

ρδc3 as a function of A/δc for monodisperse and bidisperse oblate ellipsoids. The theoretical prediction according to Eq. (6) is shown as a dotted line for comparison.

Close modal

Based on the observations above, we have analyzed some experimental data available in literature for CNT and GNP nanocomposites. As we are interested in the effect of polydispersity, we have chosen the experimental data from sources that have reported the actual size distribution of the nanofillers. We have used the data from Castillo et al.7 for multi-walled CNTs (MWCNTs) and from Tkalya et al.29 for GNPs.

The experimental data can be compared with the tunneling-percolation based theories in the following two ways:

  1. Using Eqs. (3) and (5) for polydisperse CNTs and GNPs, respectively.

  2. Using Eq. (2) based on the critical path assumption.8 

The first approach is based on continuum percolation of soft-shell hard-core fillers (or cherry-pit model), where the soft shell thickness corresponds to a tunneling-induced lengthscale λ. The second approach uses the exponential dependence of tunneling conductance (Eq. (1)) on the filler separation, but does not impart a sharp percolation cutoff and involves the tunneling lengthscale ξ, which is expected to be proportional but not necessarily equal to the lengthscale λ of the first approach.

It is clear from Figs. 2 and 9 that predictions based on the percolation connectedness theory are accurate for higher aspect ratios that are generally used in experiments. Therefore, Eqs. (3) and (5) are used to compare the model predictions with experimental data in the following. The objective is to estimate λ and ξ using the size polydispersity information along with the σϕc data from the experiments. Consistent estimations of λ and ξ within the expected range8 would indicate the importance of the polydispersity based formulation.

Following the first approach, ϕc is obtained from Eq. (3) as ϕc=D2/2λLw, where λ/2=δc/2 is tunneling soft-shell thickness. Here, Lw can be calculated from the MWCNT length distributions obtained experimentally. The polydispersity and conductivity data provided for MWCNTs from five different suppliers presented in Ref. 7 are used. We assume that all MWCNTs have the same D that is provided by the manufacturers.7 Therefore, the only unknown, λ, can be estimated using the experimentally obtained ϕc. The length distributions for the five cases are extracted from Fig. 1 of Ref. 7 and are named as A–E following the same order therein. The λ values estimated using the experimentally obtained ϕc and Lw are summarized in Table I. Interestingly, for all cases, the estimated λ falls within 7–11 nm and is comparable to the MWCNT diameter D.

TABLE I.

Estimation of the tunneling-induced lengthscale λ for cases A–E of different MWCNTs as presented by Castillo et al.7 The diameter D and percolation threshold ϕc are taken as reported in Ref. 7, while L50,L, and Lw are obtained using the MWCNT length distribution provided therein.

SamplesD (nm)L50 (nm)L (nm)Lw (nm)ϕc (vol. %)λ (nm)
7.8 735 974 1523 0.23 8.82 
10.5 770 1104 1905 0.42 6.91 
10 1341 1604 2393 0.19 10.87 
10.5 727 1027 1772 0.34 9.25 
332 382 484 0.34 10.82 
SamplesD (nm)L50 (nm)L (nm)Lw (nm)ϕc (vol. %)λ (nm)
7.8 735 974 1523 0.23 8.82 
10.5 770 1104 1905 0.42 6.91 
10 1341 1604 2393 0.19 10.87 
10.5 727 1027 1772 0.34 9.25 
332 382 484 0.34 10.82 

The effect of polydispersity is apparent when we observe the large difference between L and Lw for each case. In their study, Castillo et al.7 used the continuum percolation based formula ϕc=D/2L50 for theoretical estimation of ϕc, where L50 is the median of the MWCNT length distribution. However, such a form clearly over-predicts ϕc because the simple formula accounts neither for tunneling nor polydispersity (Fig. 10). A somewhat better match is obtained if L is used instead of L50, i.e., ϕc=D/2L which can be obtained as a special case of Eq. (3) for a monodisperse CNT distribution with lengths L, diameter D, and δc=λ=D. However, using Eq. (3) with Lw and the λ9.5 nm (obtained from the average of λ for all cases in Table I) provides excellent comparison with the experimental data, as shown in Fig. 10.

FIG. 10.

Comparison of the experimental ϕc with that of the theoretical predictions based on median L50 as ϕc=D/2L50, mean L as ϕc=D/2L, and weighted mean Lw as ϕc=D2/2λLw with λ = 9.5 nm. The unit slope line is provided to understand how well the predictions match with experimental data.

FIG. 10.

Comparison of the experimental ϕc with that of the theoretical predictions based on median L50 as ϕc=D/2L50, mean L as ϕc=D/2L, and weighted mean Lw as ϕc=D2/2λLw with λ = 9.5 nm. The unit slope line is provided to understand how well the predictions match with experimental data.

Close modal

Next, the second approach based on Eq. (2) is used to compare the model results with experimental data. For this purpose, we need conductivity σ as a function of ϕc which is provided in Fig. 2 of Ref. 7 for cases A–E. Using ϕc values, the corresponding δc is obtained using Eq. (3), i.e., δc=D2/2ϕcLw. Then, using

lnσ=2ξδc+lnσ0,
(9)

ξ can be estimated from the slope of the line fitted to points plotted on a lnσδc scale. The line is fitted only through the points above the percolation threshold where the critical path assumption is valid due to minimal effect of finite matrix conductivity.8 The points corresponding to volume fractions below the percolation threshold were observed to deviate significantly from the predicted behavior. The estimated ξ values for cases A–E (Fig. 11), summarized in Table II, are found to fall within 1.5–4 nm, which is in agreement with the expected ξ.8 Using the polydispersity information along with Eq. (3) have sufficiently narrowed down the range for ξ. The estimated ξ can be used as an input to predict σ of MWCNT nanocomposite for a given ϕc.

FIG. 11.

Estimation of tunneling lengthscale ξ through fitting using Eq. (9) for cases A–E of different MWCNTs studied in Castillo et al.7 

FIG. 11.

Estimation of tunneling lengthscale ξ through fitting using Eq. (9) for cases A–E of different MWCNTs studied in Castillo et al.7 

Close modal
TABLE II.

Estimated tunneling lengthscale ξ using Eq. (9) for cases A–E of different MWCNTs studied by Castillo et al.7 

CaseABCDE
ξ (nm) 3.2 ± 0.6 1.78 ± 0.8 3.9 ± 0.9 1.54 ± 0.9 4 ± 2 
CaseABCDE
ξ (nm) 3.2 ± 0.6 1.78 ± 0.8 3.9 ± 0.9 1.54 ± 0.9 4 ± 2 

Tkalya et al.29 used dynamic light scattering technique to obtain diameter distribution of GNPs and have shown that theoretical predictions using Eq. (5) are in qualitative agreement with experimentally obtained ϕc for polydisperse GNPs. The theoretical prediction is obtained by ignoring the thickness polydispersity and focusing only on the diameter polydispersity of the GNPs. Graphene sheet thickness L and tunneling-induced soft-shell thickness λ/2 remain as fitting parameters because both are not known a priori. However, by fitting, they have found λ/L4.5 which correctly captures the trends and qualitatively agrees with the experimental data. This is the first method of comparison as described above. Their results show that size polydispersity plays an important role in determining the percolation threshold.

Following the second method of comparison, Eq. (9) used along with Eq. (5), which relates δc with ϕc for polydisperse discs, gives

lnσ=8LD2ξϕcA+lnσ0.
(10)

Therefore, plotting lnσ as a function of 8D2/ϕcA, we can estimate ξ/L from the slope of the line fitted to the experimental data (Fig. 12) for different samples from Ref. 29. Similar to the MWCNTs case above, only points above the percolation threshold are considered to obtain the fits.8 Points below the percolation threshold can be observed to deviate from the fitted line in Fig. (12) for some cases. The diameter distributions reported therein are used to calculate A and D2. The estimated ξ/L are summarized in Table III. The GNP sheet thickness L depends on the degree of exfoliation among other factors and, although not exactly known, is expected to be between 0.3 and 2 nm.29 For the first three cases, ξ/L1, which results in an admissible estimate of ξ within 1–2 nm. For the last case, however, ξ/L=0.6, which results in a considerably smaller estimate of ξ. A similar issue with case D is reported in Ref. 29 where the experimentally obtained ϕc is much higher than the theoretical prediction from Eq. (5). For all other cases, the theoretical predictions are in excellent quantitative agreement with the experiments. Thus, by correctly accounting for the size polydispersity, accurate estimates of ξ confined to a narrower range are obtained. Using ξ values estimated above, Eq. (10) can be used to obtain conductivity of GNP nanocomposite at a given ϕc.

FIG. 12.

Estimation of ξ/L through fitting using Eq. (10) for four types of GNP dispersions as presented by Tkalya et al.29 The sample names are adopted from their article.

FIG. 12.

Estimation of ξ/L through fitting using Eq. (10) for four types of GNP dispersions as presented by Tkalya et al.29 The sample names are adopted from their article.

Close modal
TABLE III.

Estimated of ξ/L using Eq. (10) for four types of GNP dispersions as presented by Tkalya et al.29 The sample names are adopted from their article.

SamplesAA-LCA-HCB
ξ/L 0.8 ± 0.2 0.8 ± 0.4 1.1 ± 0.1 0.6 ± 0.2 
SamplesAA-LCA-HCB
ξ/L 0.8 ± 0.2 0.8 ± 0.4 1.1 ± 0.1 0.6 ± 0.2 

Both approaches used herein provide reasonable estimates of λ and ξ within the expected range. However, the tunneling lengthscale ξ in the second approach is precisely defined through Eq. (1) compared to λ in the first approach which is somewhat vague. Another estimate of λ can be obtained from the crossover point between tunneling dominated and matrix dominated regimes using fits obtained using Eq. (9). Using this apparent percolation threshold would result in lower λ values.8 However, fits based on Eq. (9) in Figs. 11 and 12 are observed to consistently over predict the effective conductivity near and below the percolation threshold. Therefore, although this apparent percolation threshold based on Eq. (9) would give lower values of λ, it may not be accurate as Eq. (9) does not seem to perform well near the percolation transition region.

It should also be noted that although polydispersity is accounted for in our model, we have neglected the effect of particle agglomerations and segregation in this study. Incorporating such effects in the model would involve quantifying the effect of anisotropic spatial distribution of fillers on the δcϕc relations. This would help in further refining the estimates of λ and ξ, and improve predictability of the models.

Theory, simulations, and experiments have all confirmed that the percolation threshold is very sensitive to the degree of polydispersity. Therefore, along with conductivity measurements, a detailed characterization of nanofiller size distribution and spatial distribution should be reported in experimental studies on conductive polymer nanocomposites.1,7,29

In summary, we have performed MC simulations of polydisperse prolate and oblate ellipsoids using the critical path based tunneling-percolation model. For polydisperse prolate ellipsoids, provided a sufficiently large L2/D, the quasi-universal behavior is obtained, as shown for polydisperse spherocylinders in Ref. 16. For polydisperse oblate ellipsoids, the data are observed to follow a quasi-universal functional form for a sufficiently large A/L. Moreover, the quasi-universal functions approach the second-virial approximation based theoretical predictions9 for rodlike and platelike fillers as the volume fraction increases.

As the connectedness percolation theory9 is in good agreement with MC results at higher aspect ratios, experimental data on polydisperse MWCNTs7 and GNPs29 are analyzed using theoretical predictions. The tunneling lengthscale ξ predicted using Eq. (2) along with the polydispersity and conductivity data from experiments is found to be well within the expected range. The estimated ξ can be used as an input to predict σ of MWCNT or GNP nanocomposites for a given ϕc.

Apart from polydispersity, other factors such as clustering of nanofillers, curling of CNTs, crumpling of GNPs, etc., are also expected to affect the percolation behavior.13 The analytical models have provided some important insights regarding the impact of such properties on percolation threshold, but are applicable, in general, for the limiting cases. Further work through modeling is in progress to understand the effect of such factors on the percolation threshold and the effective electrical conductivity of nanocomposites.

The authors gratefully acknowledge the use of the Taub cluster resources provided under the CSE program at University of Illinois. This work was financially supported by the NSF Center for Novel High Voltage/Temperature Materials and Structures [NSF I/UCRC (IIP-1362146)].

1.
R. M.
Mutiso
and
K. I.
Winey
,
Prog. Polym. Sci.
40
,
63
(
2015
).
2.
J.-M.
Thomassin
,
C.
Jrme
,
T.
Pardoen
,
C.
Bailly
,
I.
Huynen
, and
C.
Detrembleur
,
Mater. Sci. Eng., R
74
,
211
(
2013
).
3.
N.
Hu
,
Y.
Karube
,
C.
Yan
,
Z.
Masuda
, and
H.
Fukunaga
,
Acta Mater.
56
,
2929
(
2008
).
4.
Y.
Liu
,
S.
Chakrabartty
,
D. S.
Gkinosatis
,
A. K.
Mohanty
, and
N.
Lajnef
, in
IEEE Biomedical Circuits and Systems Conference, 2007, BIOCAS 2007
(
IEEE
,
2007
), pp.
119
122
.
5.
G.
Yin
,
N.
Hu
,
Y.
Karube
,
Y.
Liu
,
Y.
Li
, and
H.
Fukunaga
,
J. Compos. Mater.
45
,
1315
(
2011
).
7.
F. Y.
Castillo
,
R.
Socher
,
B.
Krause
,
R.
Headrick
,
B. P.
Grady
,
R.
Prada-Silvy
, and
P.
Ptschke
,
Polymer
52
,
3835
(
2011
).
8.
G.
Ambrosetti
,
C.
Grimaldi
,
I.
Balberg
,
T.
Maeder
,
A.
Danani
, and
P.
Ryser
,
Phys. Rev. B
81
,
155434
(
2010
).
9.
R. H.
Otten
and
P.
van der Schoot
,
J. Chem. Phys.
134
,
094902
(
2011
).
10.
V.
Ambegaokar
,
B. I.
Halperin
, and
J. S.
Langer
,
Phys. Rev. B
4
,
2612
(
1971
).
11.
W.
Bauhofer
and
J.
Kovacs
,
Compos. Sci. Technol.
69
,
1486
(
2009
).
12.
I.
Balberg
,
C. H.
Anderson
,
S.
Alexander
, and
N.
Wagner
,
Phys. Rev. B
30
,
3933
(
1984
).
13.
A. V.
Kyrylyuk
and
P.
van der Schoot
,
Proc. Natl. Acad. Sci. U.S.A.
105
,
8221
(
2008
).
14.
R. H.
Otten
and
P.
van der Schoot
,
Phys. Rev. Lett.
103
,
225704
(
2009
).
15.
A. P.
Chatterjee
,
J. Chem. Phys.
132
,
224905
(
2010
).
16.
B.
Nigro
,
C.
Grimaldi
,
P.
Ryser
,
A. P.
Chatterjee
, and
P.
van der Schoot
,
Phys. Rev. Lett.
110
,
015701
(
2013
).
17.
R. M.
Mutiso
,
M. C.
Sherrott
,
J.
Li
, and
K. I.
Winey
,
Phys. Rev. B
86
,
214306
(
2012
).
18.
Y.
Chen
,
F.
Pan
,
S.
Wang
,
B.
Liu
, and
J.
Zhang
,
Compos. Struct.
124
,
292
(
2015
).
19.
H.
Meyer
,
P. v. d.
Schoot
, and
T.
Schilling
,
J. Chem. Phys.
143
,
044901
(
2015
).
20.
J. D.
Sherwood
,
J. Phys. A: Math. Gen.
30
,
L839
(
1997
).
21.
J. W.
Perram
and
M. S.
Wertheim
,
J. Comput. Phys.
58
,
409
(
1985
).
22.
M. A.
Miller
,
J. Chem. Phys.
131
,
066101
(
2009
).
23.
M. E.
Newman
and
R. M.
Ziff
,
Phys. Rev. E
64
,
016706
(
2001
).
24.
B.
Krause
,
T.
Villmow
,
R.
Boldt
,
M.
Mende
,
G.
Petzold
, and
P.
Ptschke
,
Compos. Sci. Technol.
71
,
1145
(
2011
).
26.
S.
Wang
,
Z.
Liang
,
B.
Wang
, and
C.
Zhang
,
Nanotechnology
17
,
634
(
2006
).
27.
A. P.
Chatterjee
,
J. Chem. Phys.
141
,
034903
(
2014
).
28.
M.
Mathew
,
T.
Schilling
, and
M.
Oettel
,
Phys. Rev. E
85
,
061407
(
2012
).
29.
E.
Tkalya
,
M.
Ghislandi
,
R.
Otten
,
M.
Lotya
,
A.
Alekseev
,
P.
van der Schoot
,
J.
Coleman
,
G.
de With
, and
C.
Koning
,
ACS Appl. Mater. Interfaces
6
,
15113
(
2014
).