Monte Carlo Model for Piezoresistivity of Hybrid Nanocomposites

Large improvements in electrical conductivity have been obtained by dispersing carbon nanotubes (CNT) and graphene nanoplatelets (GNP) fillers into insulating polymer matrices [1,2]. In addition to improved electrical conductivity, filler network modification of carbon nanostructure-based nanocomposites leads to piezoresistive behavior. Several experimental studies characterize the gauge factor of CNT-based nanocomposites and evaluate their strain-sensing behavior under static and dynamic loading conditions [3–5]. On the other hand, numerical studies of this phenomenon focus on the underlying mechanisms. Behnam and Ural [6] used a quasi-three-dimensional model, to study the effect of CNT alignment, resistance ratio, length, and density on the model resistivity. Rahman and Servati modeled the impact of tunneling resistance on the nanocomposite gage factor [7]. Wang and Ye [8] considered parameters such as polymer's Poisson's ratio, CNT diameter, orientation, and density to maximize the value of average junction gap variation. Bao et al. [9,10] used an experimentally informed Weibull distribution of fillers in their network model and report that electron tunneling plays a more dominant role in the conductivity of single-wall CNT composites than using multiwall CNT fillers. Using a percolation network model, Gong and Zhu [11,12] reported that CNT deformation is responsible for the asymmetry and nonlinearity of the piezoresistivity in CNT polymer composites.

In addition to studies on CNT fillers, several recent experimental studies showed a significant improvement in electrical and mechanical properties with GNP reinforcement [13–16]. Numerical simulations of graphene-based composites pose an additional challenge because of many variations in the simulation, caused by the two-dimensional shape of the filler particles. Some previous studies dealt with percolation and conductivity with two-dimensional (2D) fillers such as circular disks and ellipses [17]. These ideas were later adapted to graphene platelets. Oskouyi et al. [18,19] used a three-dimensional Monte Carlo model combined with a finite element-based conductivity model to predict the electromechanical behavior of graphene-based nanocomposites. Yi and Tawerghi [20] evaluated the percolation threshold for 2D interpenetrating plates of circular, elliptical, square, and triangular shapes dispersed in a three-dimensional space. They observed that noncircular geometries and corner angles produce lower percolation threshold.

Given the behavior of CNT and graphene filler networks individually, a natural extension is combining these fillers to form a hybrid filler network for increased network connectivity. Hybrid nanocomposites with two or more fillers of different geometric shapes and aspect ratios have achieved better thermal, electrical, and mechanical properties in the composite materials [21–29]. An explanation of this improvement is that bridging of CNT network by planar graphitic platelets which results in a continuous network of fillers [25,27,28]. In addition, better control of dispersion can be affected by using this combination of fillers [24]. Specifically with respect to electrical conductivity, Safdari and Al-Haik [27,28] reported that low volume fraction addition of CNTs to a GNP-epoxy composites results in an order of magnitude increase in electrical conductivity and a 50% decrease in percolation threshold. Yue et al. [24] evaluated CNT/GNP mix ratios that improve conductivity and reduce the percolation threshold.

Despite the numerous studies focused on the use of CNTs or GNPs as fillers to improve the electrical conductivity and electromechanical behaviors of polymer composites, relatively few studies have examined hybrid CNT–GNP composites. In particular, piezoresistive behavior of CNT–GNP hybrid nanocomposites has not been numerically investigated, to the best of our knowledge. In the present study, we have developed a two-dimensional Monte Carlo model to estimate the electrical conductivity and the piezoresistivity of nanocomposites, with microstructures consisting of elliptically shaped GNPs, mixed in with a two-dimensional CNT network. Stochasticity is naturally inherent in these microstructures; consequently, there is uncertainty in the constituent model parameters. We study the uncertainty in the model outputs by simulating a large number of sample paths (Monte Carlo simulations) and address the uncertainty due to model assumptions as well as parametric uncertainty due to the input parameters in simulating the microstructures. Some of the parameters considered include filler geometry, aspect ratio, orientation, CNTs to GNPs volume fraction ratio, and mechanical loading. The effect of these parameters on the percolation threshold, electrical conductivity, and piezoresistivity is examined with the objective of designing responsive microstructures of hybrid nanocomposites.

where $li$ and $θi$ are the length and polar angle of the ith CNT, respectively, as shown in Fig. 1(a). $x1i$⁠, $y1i$⁠, and $θi$ follow a uniform distribution, while $li$ follows a Weibull distribution such that

where “rand” indicates uniformly distributed random numbers in the interval $]0,1[$ and $(a,b)$ are the parameters of Weibull distribution. The Weibull parameters are determined by the nominal and standard deviation of experimental CNT lengths [30]. Some of the earlier modeling studies consider CNTs of uniform length [7,8,31], however, including the length distribution according to experimental parameters results in a more realistic representation of the microstructure.

This CNT generation procedure was performed continuously until the desired CNT volume fraction was reached. Even though CNTs are represented as one-dimensional structures, a fixed diameter D is assigned to the CNTs to facilitate computing the intersections and volume fraction. Here, the volume fraction is simply defined as the ratio of the total area of all the CNTs (⁠$li×D)$ and the area of RVE (⁠$Lx×Ly)$⁠. Periodic boundary conditions (PBC) [31] were used to reposition the CNTs when the endpoints of CNTs generated by Eq. (1) were located outside of the RVE.

Next, we consider the incorporation of elliptical graphene platelets to form a hybrid filler network. Most studies on percolation of graphene-based nanocomposites modeled graphene as circular particles because of the simplicity in geometry definition, and ease in computation of the distance between particles [18,19,27]. Note that the two fillers generated in this work are in plane with the 2D RVE.

The elliptical particles were added to the RVE until the desired graphene volume fraction was reached. Equation (3) suggests that the center of the ellipse is always located inside the RVE. However, the ellipse may intersect with one or two of the RVE boundaries. Periodic boundary conditions were used in such cases. Note that the implementation of PBCs for elliptical particles has additional steps compared to circular particles [32]. If an ellipse intersects with a single boundary segment, one additional ellipse $(x0i+1,y0i+1,θ0i)$ was compensated (Fig. 2(a)). If the ellipse intersects with two boundary lines, on or outside the segments, two additional ellipses $(x0i+1,y0i+1,θ0i)$ and $(x0i+2,y0i+2,θ0i)$ were compensated (Figs. 2(b) and 2(c)). If one of the RVE corner is located inside the ellipse, three additional ellipses $(x0i+1,y0i+1,θ0i)$⁠, $(x0i+2,y0i+2,θ0i)$⁠, and $(x0i+3,y0i+3,θ0i)$ were compensated (Fig. 2(d)). Figures 2(b) and 2(c) can happen only with radially asymmetric shapes like ellipses. As a modeling constraint, the ellipses are impenetrable, i.e., not allowed to overlap. This constraint is needed to make the 2D model closer to a three-dimensional microstructure.

In order to reduce computational effort, we used an algorithm wherein a small number n_p of equidistant points (⁠$xkj, ykj$⁠) on ellipse j satisfy Eq. (7), with $k=1,2,…,np$⁠.

Percolation network conductivity models use either soft core or hard core approaches to determine the junctions between filler particles. Soft core approach assumes fully permeable fillers wherein geometric and electrical contact occur simultaneously [31]. In hard core approaches, a soft shell encompasses impenetrable hard core, and tunneling intersections occur in the soft shell. Most of the 2D numerical studies [6,7] as well as many three-dimensional modeling studies [9–12], use a soft-core model for determining junctions between fillers. Berhan and Sastry [33] showed that soft core model reduces the contribution of tunneling, due to the overestimation of the number of chemically bonded junctions and thereby introduces non-negligible error in percolation modeling. Hard core approach is more accurate and realistic but also computationally more intensive, especially for fillers with complex geometries.

In our model, we used a modified soft core approach for CNT–CNT and CNT–GNP interactions. Tunneling effect happens between two filler particles when the shortest distance between them is less than the maximum effective distance (⁠$dcutoff$⁠) for tunneling effect (see Fig. 3(b)). The geometrical contact between CNTs in the model was transformed into a tunneling junction for resistance calculation. For any CNT–CNT intersection, a random tunneling distance, $dmn$ such that $(D≤dmn≤D+dcutoff)$ was generated, to create two distinct contact nodes $m$ and $n$ (one on each of the two CNTs in contact) such that the distance between those points (⁠$m$ and n) was equal to the generated tunneling distance $dmn$ (see Fig. 3(a)). This tunneling distance is used for percolation and resistance calculations, thereby reducing the error in softcore method due to underestimation of tunneling junctions, without adding significant computation. A similar approach was used for the intersection between GNPs and CNTs as shown in Fig. 3(e).

The parameters of ith ellipse were obtained by replacing $a0$ and $b0$ in Eq. (6) by $aei$ and $bei$⁠. Then, overlapping between the ith ellipse's tunneling softcore and the jth ellipse was checked using $np$ points (⁠$xkj, ykj$⁠) and applying Eq. (7), with $k=1,2,…,np$⁠. In the case of overlapping, the shortest distance between (⁠$xkj, ykj$⁠) points on the jth ellipse and (⁠$xki, yki$⁠) points on the ith ellipse (not its soft core) were used to generate a tunneling resistor segment with a tunneling distance $dmn$⁠, as shown in Fig. 3(d). The tunneling segments between the filler particles along with the filler particles were used to create clusters that lead to percolation networks. Note that all the nodes created due to tunneling on GNPs are projected on the GNP resistor segment (see Figs. 3(d) and 3(e), from left to right).

where $ljk$ is the CNT length between the contact points $j$ and $k$⁠, $D$ is the diameter of the CNT, and $σcnt$ is the CNT intrinsic electrical conductivity. In this study, we used an experimentally determined effective mean conductivity as $σcnt$ [11]. Note that the electrical conductivity of the polymer matrix is neglected.

where $me$ is the mass of electron, $ΔE$ (1 eV) is the height of the barrier (the difference of the work functions between the CNT and the polymer), $dvdw$ is the van der Waals separation distance. Examples of resistance network between fillers are shown in Figs. 3(c)–3(e), with the different type of electrical resistances.

where $np$ is the number of simulations with the existence of at least one conductive path, in a total of $N$ simulations. The electrical percolation threshold corresponds to a percolation probability of $0.5$⁠.

Once a conductive network, which spans both electrodes was found, the segments that did not participate in conducting current were removed for computational efficiency using the Dulmage–Mendelsohn decomposition method [35]. Then, the remaining conductive network was transformed into a resistor network by calculating the resistance between the different nodes (⁠$RintCNT=((4ljk)/(σcntπD2))$⁠, $RintGNP=(R/L)Lik$⁠, and $Rc=(h/(2e2MT))$⁠). The final step was to calculate the resistance of that percolating network from the positive semidefinite matrix equations representing Kirchhoff's current and Ohm's laws [9,36]. Cholesky decomposition algorithm for sparse matrices [36] was used to solve these matrix equations and find the conductance of the RVE. Monte Carlo simulations on several representative rectangles were performed to obtain the averaged conductance of the networks.

Several experimental studies confirm that nanocomposites with CNTs or graphene particles exhibit a change in resistance when subject to mechanical deformation [3–5]. We study this phenomenon by subjecting the RVE to uniaxial tensile deformation. Consider a CNT–graphene nanocomposite under an incremental uniaxial strain $Δε$ along the y-axis of the RVE. Any applied strain resulted in a change in position and orientation of the particles. The re-orientation model for the piezoresistivity for CNTs in Ref. [37] was applied for both CNTs and graphene particles in this study. In the re-orientation model described in Eqs. (20) and (21), a uniform strain field and a constant Poisson's ratio is assumed in the polymer. Also a perfect interface between the polymer matrix and the fillers is assumed. Hence, all points in the RVE including fillers are assumed to deform under mechanical strain with elastic properties of the matrix. Since the intrinsic piezoresistivity of the fillers are not considered in our calculations, the effects of this assumption on the accuracy especially for small strains are expected to be very small.

We use Monte Carlo simulations based on the model described above to compute the percolation probability and the electrical conductivity, first with only CNTs as fillers and then with CNTs and GNPs as fillers. The different parameters (size, content, and aspect ratio) of GNP as well as CNT content are kept constant to specific set of values, while the position and alignment of GNP and CNT as well as the length of CNT are randomly generated using Monte Carlo simulations, as shown in Eqs. (2) and (3).

To have a conductive nanocomposite, after the application of strain, the percolation network should remain after deformation. Our simulations show that when the CNT volume fraction is less than 0.10, the percolation network becomes discontinuous (i.e., the percolation path does not hold all the time) at a strain of 0.6% or less. Therefore, to study the piezoresistivity of hybrid composites, the CNT volume fraction was fixed at 0.10 in all the piezoresistivity simulations, and only the GNP parameters were varied for determining optimal microstructures. This ensures that a percolated network of CNTs per se exists before and after the application of strain.

The simulation parameters for the fillers are tabulated in Table 1. Statistical variation in microstructures is reduced by averaging over a large number of randomized microstructures. RVE cell of size 25 μm × 25 μm is used. Figure 4 shows the convergence study over 10,000 Monte Carlo simulations. We can observe convergence of the electrical conductivity after 2000 simulations. A minimum of 2000 Monte Carlo simulations were performed for each variation of geometric features or fillers loading for both conductivity and piezoresistivity calculation.

In order to investigate the effect of GNP content on the behavior of hybrid nanocomposites, the GNP size was kept constant, while the content of GNP (controlled by the GNP-to-CNT volume fraction ratio) was varied. GNP-to-CNT volume fraction ratios (GNP/CNT) of 0.5, 1, and 2 were considered. Three different elliptical graphene dimensions were studied: the ellipses with semimajor and minor axes, a = 4 μm and b = 2 μm, a = 2 μm and b = 1 μm, and a = 1 μm and b = 0.5 μm named, respectively, as GNP (4,2), GNP (2,1), and GNP (1,0.5). Figure 5(a) reveals that the addition of GNP (4, 2) to the CNT nanocomposites, for any GNP-to-CNT volume fraction ratio, increased the percolation threshold of the hybrid nanocomposites, compared to nanocomposites with only CNT fillers. Moreover, an increase of the GNP-to-CNT volume fraction ratio led to a steady increase of the percolation threshold.

Figure 5(b) shows that the addition of GNP (4, 2) for any GNP content increased the conductivity of the nanocomposites. It also shows that a higher GNP-to-CNT volume fraction ratio led to a higher conductivity of the nanocomposites. Since the size of the GNP was constant while the GNP-to-CNT ratio was increased, more GNP particles were being added to the nanocomposites. Increase of the number of GNP particles might have increased the possibility of forming more junctions between the fillers.

Figure 5(c) shows that the addition of GNPs to the nanocomposites increased the piezoresistive performance of the nanocomposites for all three graphene contents, compared to the CNT monofiller nanocomposites. The simulations for the piezoresistivity were all performed for a CNT volume fraction of 0.10. Figure 5(c) also shows that when the content of GNP was increased the piezoresistivity also increased for all the values of applied mechanical strains. It means that the hybrid nanocomposite with increased GNP content is more sensitive to the applied strain. A similar trend for percolation, conductivity, and piezoresistivity was observed for the other two graphene sizes (GNP (2, 1) and GNP (1, 0.5)). It is known that the piezoresistivity of CNT-based nanocomposites with lower CNT content is more sensitive to the applied strain [11]. This limits the gage factor because of the limit imposed by the percolation threshold. Our result indicates that the gage factor of hybrid CNT–GNP nanocomposites does not have that limitation (the GNP content can be increased for better piezoresistivity). Also the addition of GNP (4, 2) at a GNP/CNT ratio of 2, doubled the piezoresistivity of the hybrid nanocomposite compared to a pure CNT nanocomposite.

This section investigates the impact of the size of GNP particles on the hybrid nanocomposite. The GNP aspect ratio and content were kept constant. Only the size (surface area) of the GNP particles was changed.

The three sizes of graphene particles mentioned previously (range of semimajor axis from 1 μm to 4 μm), were used in the simulations. The electrical conductivity of the nanocomposites was determined for varying CNT volume fractions. The results in Fig. 6(a) show that an increase of the size of GNP led to a steady decrease of the percolation threshold of the hybrid nanocomposites. Also, Fig. 6(b) reveals that nanocomposites with bigger graphene particles have higher conductivity. Bigger GNP fillers in smaller number worked better than smaller particles in larger number. Hence, the effect of GNP size on the conductivity of hybrid nanocomposites is similar to the effect of CNT length on the conductivity of CNT-based nanocomposites [11].

The probability of forming more junctions between the fillers seems to increase with larger particles. The piezoresistivity behavior of the hybrid nanocomposites was modeled for a fixed CNT volume fraction of 0.10, in all the simulations. Figure 6(c) shows that the increase of the GNP particles size increased the piezoresistivity of the hybrid nanocomposites. A similar trend for percolation, conductivity, and piezoresistivity was seen for the two other GNP-to-CNT volume fraction ratio values (0.5 and 1).

To investigate the impact of the aspect ratio of the GNP particles on the hybrid nanocomposites, the GNP particles size (surface area) and content were kept constant, while the aspect ratio of the GNP was varied. The simulations were performed for two aspect ratios (2 and 8) with GNP (2, 1) and GNP (4, 0.5).

The percolation threshold decreased (Fig. 7(a)) with an increase of the electrical conductivity (Fig. 7(b)) when the GNP aspect ratio of the hybrid nanocomposites is increased from 2 to 8. The simulations on the piezoresistivity of the hybrid nanocomposites, in Fig. 7(c), with a CNT volume fraction of 0.10, show that the piezoresistivity increased when the aspect ratio of GNP was increased. Since the size of the GNPs was the same, we have the same number of GNPs in the RVE for the two aspect ratios.

We now analyze the effect of aspect ratio for varying GNP size. The major axis length and the content of the GNP particles were kept constant for two aspect ratios, 2 and 8 corresponding to GNP (4, 2) and GNP (4, 0.5) fillers. Figure 8 confirms that a higher aspect ratio leads to a better percolation, electrical conductivity, and piezoresistive performances. However, since the major axis length of the two type of GNP were all equal, we see that the GNP particle with a bigger size (GNP (4,2)) did not lead in this case to the higher conductivity or piezoresistivity (unlike in Sec. 3.2, where the aspect ratio was constant).

The electromechanical properties appear to exhibit a higher sensitivity to aspect ratio compared to GNP size. This can be a useful pointer in designing microstructures for optimal behavior. The trend in Figs. 7 and 8 was observed for the two other GNP-to-CNT volume fraction ratio values (0.5 and 1).

It has been shown in previous studies that CNT alignment could have a huge effect on the piezoresistivity of CNT-based nanocomposites [39]. In this section, we investigate the influence of the alignment of GNP particles on the piezoresistivity of the hybrid nanocomposites. First, the GNP (4, 0.5) fillers were randomly generated with different orientations with the angle bounded by $θmax$⁠, i.e., graphene particles were oriented at an angle between 0 and $θmax$⁠. Figure 9(a) shows the effect of the alignment of GNP (4, 0.5) on the piezoresistivity of the hybrid nanocomposites for values of $θmax$ (30 deg, 60 deg, 90 deg, 120 deg, and 180 deg). It also shows the direction of measurement of the conductivity and of the fillers orientation. For small values of $θmax$⁠, we observed lower piezoresistivity. The resistance ratio reached its minimum for values of $θmax$ of 30 deg and 60 deg. There is an increase in piezoresistivity at higher angles of $θmax$⁠, but there is no linear or monotonic correlation between piezoresistivity and $θmax$⁠. Note that smaller angle results in more transversely oriented particles, while larger angles denote a more random orientation. We however observe a significant increase (five times compared to $θmax$ of 30 deg) in piezoresistivity when all GNPs were oriented longitudinally, i.e., at a 90 deg angle.

A second set of simulations for GNP alignment was performed where GNP (4, 0.5) were randomly generated with orientations chosen from equal adjacent ranges of angle, of amplitude 30 deg as shown in Fig. 9(b). We observed that angle ranges adjacent to 90 deg gives the highest piezoresistivity values (60–90 deg and 90–120 deg) while the angle ranges adjacent to 0 deg or 180 deg gives the lowest piezoresistivity values (0–30 deg and 150–180 deg). Moreover, the piezoresistive performance of the hybrid nanocomposite in function of GNP orientation seems to be symmetric to 90 deg orientation and decreases monotonically when we move away from 90 deg. We find a similar trend when the simulations were repeated with GNP (4, 2) fillers.

Figure 10 shows the different parameters of GNP examined as well as their effect on the piezoresistivity of the hybrid nanocomposites with CNT and GNP content kept the same for all microstructures. The resistance change ratio of a nanocomposite purely based on CNT at 0.6% of strain is shown to be about 0.32%, with a corresponding gage factor of 0.5. With the addition of the smallest size of graphene particle examined in this study, GNP (1, 0.5) to the CNT nanocomposites, the gage factor of the hybrid nanocomposites increased 1.6 times. Gong and Zhu, showed that the addition of more CNT (instead of graphene particles) to the initial CNT nanocomposites decreased the gage factor of the nanocomposite, because an increase in CNT volume fraction leads to the decrease of the rate of change of junction numbers between fillers after strain was applied [11]. Thus, the improvement of the piezoresistivity of the hybrid nanocomposite means that the rate of change of the junction number between fillers after strain, might have increased, with the addition of GNP particles. This might be due to the fact that GNP as two-dimensional particles are able to make more junctions with neighboring fillers in the RVE, compared to one-dimensional CNTs.

The next microstructure simulated was with bigger graphene particles (GNP (4, 2)), with the same graphene content as before. A gage factor 1.5 times greater than with the smaller GNP was observed. This improvement with a bigger particle size might be explained by the fact that, with a network of CNTs already existing, smaller particles, even though in bigger number are more likely to fall in between the CNTs or to overlap CNTs that do not participate in the percolated network. However, bigger GNP particles, even in smaller number, are more likely to overlap more neighboring fillers, some of which will be part of the percolated network. This allows the bigger GNP particles to increase, more effectively the number of fillers in the percolated network, compared to smaller GNP particles.

A microstructure with the same content of GNP, but with higher aspect ratio (GNP (4, 0.5) fillers, is compared to the previous microstructure. An increase of the gage factor up to 1.7 times was observed for the same strain of 0.6%. This indicates that piezoresistivity is more sensitive to aspect ratio than to size of the particles. Finally, GNP (4, 0.5) fillers, at the same GNP content as in the previous microstructure are used but all oriented parallel to the loading direction. The gage factor is shown to double when compared to the randomly oriented microstructure of GNP (4, 0.5). This might be due to the fact that with an existing percolated network of CNTs, GNP fillers oriented in the direction of the conductivity measurements are able to form junctions and span the existing CNT fillers to a greater extent.

Effect on the electrical conductivity of the same microstructures is shown in Fig. 11. In these simulations, CNT volume fraction is also varied along with GNP parameters. Figure 11 shows that the improvement of the percolation threshold, and of the electrical conductivity, when the GNP parameters are varied, follows the same trend as the one of the piezoresistivity. The reason behind those improvements are the same as the one mentioned for the piezoresistivity, i.e., rate of change of the number of CNT–GNP tunneling junctions.

Ranking the parameters, we find that GNP alignment (all 90 deg orientation) had the highest influence resulting in 5.94 times increase in gage factor compared to CNT monofiller composite. Other factors in the decreasing order are aspect ratio, GNP size, and simply the effect of adding the second phase (GNP) enhanced the gage factor up to 3.75, 2.21, and 1.83 times that of CNT monofiller composite, respectively.

We developed a new CNT–GNP percolation network model, and used it to investigate the mechanisms of electrical percolation, conductivity, and piezoresistivity of hybrid CNT–GNP nanocomposites, where tunneling is the electron transport mechanism. This model considers for the first time, the effect of addition of GNP on the piezoresistivity behavior of hybrid nanocomposites. The GNP fillers are modeled as ellipses in the Monte Carlo simulations. We find that the addition of GNP particles to CNT nanocomposites enhances significantly the electrical conductivity and piezoresistive behavior of the nanocomposites. Furthermore, parametric analysis has been conducted to investigate the impact of different GNP parameters on the piezoresistivity of the hybrid nanocomposites. We show that: (i) a conductive network with higher GNP loadings, has a higher percolation threshold, and a higher electrical conductivity and piezoresistivity, (ii) an increase of the GNP size or aspect ratio decreases the percolation threshold and increases the conductivity and piezoresistivity of the hybrid nanocomposites, (iv) alignment of the GNP particles has a big effect on both the piezoresistivity and electrical conductivity, and (v) GNPs uniformly aligned in the direction of electrical conductivity measurements (or in a direction close to that) lead to a significantly higher piezoresistive behavior, up to six times greater than that of nanocomposites based on only CNT. We also attempt to elaborate on the reasons behind the improvement of piezoresistivity behavior observed in the different microstructures in this study.

• ERAU seed grant.

• NASA SBIR grant.