Analytical modeling of the electrical conductivity of CNT-filled polymer nanocomposites

Abstract

Electrical conductivity of most polymeric insulators can be drastically enhanced by introducing a small volume fraction $(~1\%)$ of conductive nanofillers. These nanocomposites find wide-ranging engineering applications from cellular metamaterials to strain sensors. In this work, we present a mathematical model to predict the effective electrical conductivity of carbon nanotubes (CNTs)/polymer nanocomposites accounting for the conductivity, dimensions, volume fraction, and alignment of the CNTs. Eshelby’s classical equivalent inclusion method (EIM) is generalized to account for electron-hopping—a key mechanism of electron transport across CNTs, and is validated with experimental data. Two measurements, namely, the limit angle of filler orientation and the probability distribution function, are used to control the alignment of CNTs within the composites. Our simulations show that decreasing the angle from a uniformly random distribution to a fully aligned state significantly reduces the transverse electrical conductivity, while the longitudinal conductivity shows less sensitivity to angle variation. Moreover, it is observed that distributing CNTs with non-uniform probability distribution functions results in an increase in longitudinal conductivity and a decrease in transverse conductivity, with these differences becoming more pronounced as the volume fraction of CNTs is increased. A reduction in CNT length decreases the effective electrical conductivity due to the reduced number of available conductive pathways while reducing CNT diameter increases the conductivity.

1. Introduction

The discovery of carbon nanotubes (CNTs) in 1991 by Iijima [1] was an influential achievement because of their strength, low weight, and high mechanical, thermal, and electrical properties [2, 3]. To take advantage of these extraordinary physical properties, CNT-reinforced composites are developed as an emerging class of advanced lightweight nanocomposites [4–7]. The properties and behavior of various structures made of these composite materials have been widely evaluated and reported in the literature using different experimental, analytical, and numerical approaches, and numerous potential applications spanning across automotive, electronics, energy storage, biomedical engineering, sports equipment, and structural reinforcement are reported for this class of composites [8–14].

Analytical micromechanics theories such as the rule of mixture [15–17], Eshelby’s equivalent inclusion method (EIM) [18–20], Halpin–Tsai method [21–24], Lewis–Nielsen method [25, 26], and Mori–Tanaka method [27–29] have been widely used in many studies to evaluate the overall physical properties of composite materials with reasonable accuracy. Developed by Eshelby [20], one of the most popular theories, the EIM shows that the elastic fields inside an ellipsoidal inhomogeneity can be assumed to be uniform if the ellipsoidal particle is perfectly bonded to an infinitely extended matrix with a uniform load applied at infinity. One of the main advantages of EIM is that the solution is limited to a system of algebraic equations, and it can be applied to the different behavior of materials such as elastic-plastic, viscoelastic, and creep [30]. The EIM was extended by Dunn and Taya [31] for the problems with the coupled electro-elastic behavior of piezoelectric composites based on the rigorous electro-elastic solution of an ellipsoidal inclusion in an infinite piezoelectric medium. The method in steady-state heat conduction was extended further for the randomly oriented particle composites with uncoupled thermal and electromagnetic behavior by Hatta and Taya [32, 33]. They considered two distribution functions for the inclusions and showed that their model is able to consider the effects of interaction between different inclusion orientations. Chen and Wang [34] proposed an analytical model based on the Mori–Tanaka mean field theory and the EIM to evaluate the thermal conductivity of the composite materials. They developed a new distribution function called Kacir’s single-parameter exponential function to model the inclusion orientation. Odegard and Gates [35] developed a method to modify the constitutive models of polymer composite reinforced by randomly oriented CNTs. They captured the effect of the discrete nature of the atomic interactions at the nanoscale and the interphase between CNTs and the polymer matrix using the EIM. Three different patterns for the dispersion of CNTs including CNTs with an axisymmetric orientation, aligned, and randomly oriented CNTs were considered in their work. Seidel and Lagoudas [36] used Mori–Tanaka, self-consistent, and composite cylinders micromechanical models in conjunction with Eshelby’s method to analyze the elastic behavior of nanocomposites reinforced by CNTs. They used a tessellation procedure and a multi-phase Mori–Tanaka approach to capture the effect of clustering of CNTs in the polymer matrix. Moreover, the effects of interphase regions were taken into account. Jin et al. [37] proposed a closed-form solution for the stress caused by the disturbance of an elliptical inhomogeneity in an infinite isotropic elastic plane. Their study showed that the EIM is a useful tool at least for solving the classical elliptical inhomogeneity problem.

Several analytical models were developed to analyze the electrical conductivity and piezoresistivity of CNT-reinforced composites. In another work by Seidel and Lagoudas [38], a micromechanical model based on the Mori–Tanaka approach and EIM was proposed to investigate the effects of electron hopping by considering conductive interphase layers around particles. Their results showed that the thickness of the electrical tunneling interphase layer relative to the CNT radius gives a distinct percolation concentration in which the well-dispersed CNT fillers are in close vicinity and electrical tunneling easily happens. A simple analytical model was presented by Deng and Zheng [39] to predict the electrical conductivity of CNT-reinforced composites considering the effects of the percolation, conductive networks, conductivity anisotropy, and waviness of CNT inclusions. Their results revealed that the waviness of CNTs has an important influence in evaluating electrical conductivity. Following their model, Takeda et al. [40] proposed an analytical model for evaluating the electrical conductivity of CNT-reinforced polymers by taking into account the nanoscale effects and the electron hopping. They also conducted experiments and compared the results with their analytical model. Feng and Jiang [41] developed a hybrid analytical model to evaluate the electrical properties of CNT-reinforced polymers by incorporating electron hopping and conductive networks. An interphase layer surrounding the CNT was used to capture the nanoscale effect of electrical tunneling based on the electron hopping theory. They illustrated that electron hopping and conductive networks contribute to the electrical conductivity, while the conductive networks’ effect is dominant above the percolation threshold. Garcia-Macias et al. [12] developed a micromechanics model which can take into account the non-straightness by a helical waviness model and non-uniform dispersion of CNTs by a two-parameter agglomeration approach for evaluating the effective electrical conductivity of cement-based composites. The electron hopping and conductive networks were modeled via an interphase surrounding the CNTs and changes in the CNTs aspect ratios, respectively. Haghgoo et al. [42] studied the effects of aspect ratio, volume fraction, clustering, and distribution of CNTs on the electrical conductivity of hybrid polymer composites reinforced by CNTs and carbon fibers. The effect of electron tunneling was considered by modeling an interphase layer. They showed that the electrical conductivity highly depends on the aspect ratio, alignment angle, and the distance among the CNT particles. A micromechanical model based on the Eshelby–Mori–Tanaka approach was conducted by Mora et al. [43] to evaluate the electrical conductivity of polymer nanocomposites with agglomeration and segregation of CNT particles. Their model is able to predict the effect of segregation on the electrical conductivity of nanocomposites. Their model was validated by comparing the predicted electrical conductivity to the experimental values. Chanda et al. [44] developed two analytical models based on micromechanics and a simplified fiber–fiber contact network to calculate the electrical conductivity of random and aligned composites. They investigated the effect of aspect ratio, waviness, orientation, interphase, and percolation threshold on the effective conductivity of composites. Shingare and Naskar [45] investigated the micromechanical modeling of fiber-based composites, emphasizing the impact of random fiber orientation on piezoelastic properties. They estimated the effective properties of composites considering the different orientations, shapes, and geometrical aspects of graphene fibers. Tang et al. [46] presented a simple analytical model to estimate the percolation threshold and electrical conductivity of CNT-reinforced composites by considering the effects of waviness, dispersion, volume fraction, and size of particles. They validated their model by comparing their results with experimental values in the literature. Saberi et al. [47] presented a multi-step analytical model to predict the electrical conductivity of the short carbon fiber-reinforced hybrid composites with graphene nanoplatelet particles. They investigated the effects of volume fraction, interphase characteristics, and fiber aspect ratio. Haghgoo et al. [48] presented a two-step analytical approach considering percolation network and electron tunneling to predict the electrical resistivity and percolation threshold in CNT-reinforced hybrid nanocomposites with carbon black. Quinteros et al. [49] proposed a computational framework for evaluating the electromechanical behavior of CNT-reinforced composites during fracture by combining electrical-deformation-fracture finite element (FE) modeling with a mixed micromechanics formulation. In their model, the electrical conductivity has a nonlinear relationship with the volume content of distributed CNTs.

To summarize the above discussion, the idea of equivalent inclusion proposed by Eshelby [20] for elasticity has been widely extended to evaluate the other overall physical properties of inclusion–matrix materials, e.g., thermal conductivity [32–34, 50] and electrical conductivity [12, 38, 41, 43] of composites. However, in the case of electrical conductivity, the lack of detailed formulation derivation is notable. Calculations based on the EIM are highly sensitive to input parameters which are often overlooked. Moreover, while most analytical calculations focus on isotropic composites, CNTs are commonly distributed in a non-uniform manner in composites, influenced by factors such as the manufacturing process or deformation due to large strain. Specific applications may call for deliberate manipulation of inclusion distribution to enhance material performance and achieve optimal designs. To address these issues, we present a rigorous analytical formulation for computing effective electrical conductivity of polymer nanocomposites accounting for both conductive networks and electron hopping [32, 38, 51]. Our calculations demonstrate a high sensitivity of results on input parameters such as energy barrier for electron hopping, intrinsic CNT conductivity, and percolation threshold. Often, the exact values of these parameters are not reported in experiments and care must be taken to use appropriate values in computational models. Finally, we model the non-uniform distribution of CNTs using a limit angle of filler orientation and a probability distribution function. The EIM of homogenization is reviewed in section 2 and is generalized to include electron hopping in section 3. The model is validated with experimental data in section 4. We also discuss the role of CNT volume fraction, orientation, distribution, and dimensions in this section with the help of numerical examples. The conclusions are reported in section 5.

2. EIM

This section summarizes the EIM for calculating the effective electrical or thermal conductivity of two-phase composites [20, 32, 33, 52, 53]. The inclusions can be modeled with or without an interphase coating with a different conductivity [41, 54]. First, the method is applied to a matrix with a single ellipsoidal particle and then extended to a matrix with numerous ellipsoidal particles, considering the interactions between them.

2.1. Overall electrical conductivity of two-phase composites

In general, the overall electrical conductivity of any two-phase composite, $K_{\mathrm{ij}}$ , can be estimated by average electric current density $\overline{J_i}$ , and average electric field of composite $\overline{E_i}$ , using Ohm’s law as:

$$\overline{J_i}=K_{\mathrm{ij}}\overline{E_j},$$

(1)

where horizontal line denotes the averaged quantity over the volume $V$ as $\overline{\square }=\frac{1}{V}\int \square \mathrm{d}V$ . Unless otherwise stated, index notation is used through this section defining Einstein’s summation convention as $a_i{b}_i=\sum _{i=1}^{\dim }{a}_i{b}_i$ , where $\dim$ denotes the dimension of the geometry. The whole domain of composite represented by $\mathfrak{D}$ , can be divided by the domain of matrix ${\mathfrak{D}}^m$ and the domain of particles ${\mathfrak{D}}^p$ such that ${\mathfrak{D}}^m\cup {\mathfrak{D}}^p=\mathfrak{D}$ , and ${\mathfrak{D}}^m\cap {\mathfrak{D}}^p=\varnothing$ . By decomposing the average quantities $\overline{J_i}$ and $\overline{E_i}$ into a matrix phase and a particle phase, they can be written as:

$$\overline{J_i}=\frac{1}{V_{\mathfrak{D}}}[{\int }_{{\mathfrak{D}}^m}{J}_i\mathrm{d}V+{\int }_{{\mathfrak{D}}^p}{J}_i\mathrm{d}V],$$

(2a)

$$\overline{E_i}=\frac{1}{V_{\mathfrak{D}}}[{\int }_{{\mathfrak{D}}^m}{E}_i\mathrm{d}V+{\int }_{{\mathfrak{D}}^p}{E}_i\mathrm{d}V].$$

(2b)

Ohm’s law in each phase is given as:

$$J_i=K^m{\delta }_{\mathrm{ij}}{E}_j\mathrm{in}{\mathfrak{D}}^m,J_i=K_{\mathrm{ij}}^p{E}_j\mathrm{in}{\mathfrak{D}}^p,$$

(3)

where $K^m$ is the electrical conductivity of the isotropic matrix and $K_{\mathrm{ij}}^p=\mathrm{diag}(K_{11}^p,K_{22}^p,K_{33}^p)$ is the anisotropic electrical conductivity of particle, and ${\delta }_{\mathrm{ij}}$ is the Kronecker delta. Equation (3) can be integrated over the respective domains to give:

$${\int }_{{\mathfrak{D}}^m}{J}_i\mathrm{d}V=K^m{\delta }_{\mathrm{ij}}{\int }_{{\mathfrak{D}}^m}{E}_j\mathrm{d}V,$$

(4a)

$${\int }_{{\mathfrak{D}}^p}{J}_i\mathrm{d}V=K_{\mathrm{ij}}^p{\int }_{{\mathfrak{D}}^p}{E}_j\mathrm{d}V.$$

(4b)

By calculating $\overline{J_i}$ from equations (2a), (4a), and (4b), equation (1) yields:

$$K_{\mathrm{ij}}\overline{E_j}=\frac{1}{V_{\mathfrak{D}}}[K^m{\int }_{{\mathfrak{D}}^m}{E}_i\mathrm{d}V+K_{\mathrm{ij}}^p{\int }_{{\mathfrak{D}}^p}{E}_j\mathrm{d}V].$$

(5)

By substituting ${\int }_{{\mathfrak{D}}^m}{E}_i\mathrm{d}V$ from equation (2b), the above equation gives the overall electrical conductivity as:

$$K_{\mathrm{ij}}\overline{E_j}=K^m\overline{E_i}+\frac{1}{V_{\mathfrak{D}}}[K_{\mathrm{ij}}^p-K^m{\delta }_{\mathrm{ij}}]{\int }_{{\mathfrak{D}}^p}{E}_j\mathrm{d}V.$$

(6)

Hence, in order to estimate $K_{\mathrm{ij}}$ , we need to evaluate $\overline{E_i}$ and ${\int }_{{\mathfrak{D}}^p}{E}_i\mathrm{d}V$ in equation (6).

2.2. Matrix with a single ellipsoidal particle

Consider a single ellipsoidal particle embedded in an infinite matrix while the constant electric current density $J_i^0$ is applied at the far field as shown in Figure 1(a). By decomposing the electric field in the composite as $E_i=E_i^0+E_i^d$ , Ohm’s law inside the particle domain holds:

$$J_i=K_{\mathrm{ij}}^p{E}_j=K_{\mathrm{ij}}^p[E_j^0+E_j^d]\mathrm{in}{\mathfrak{D}}^p,$$

(7)

where $E_i^d$ is the disturbed electric field due to the existence of inhomogeneity ${\mathfrak{D}}^p$ , and $E_i^0$ is the uniform electric field due to the current density $J_i^0$ in the absence of the particle, i.e.,

$$J_i^0=K^m{E}_i^0\mathrm{in}\mathfrak{D}.$$

(8)

Figure 1.

(a) A single ellipsoidal particle embedded in an infinite matrix subjected to the applied current density and (b) the equivalent inclusion inducing transformation electric field.

As shown in Figure 1(b), consider an imaginary subdomain ${\mathfrak{D}}^p$ called an inclusion which undergoes transformation electric field but its electrical conductivity is equal to the electrical conductivity of the matrix. It can be shown that the EIM captures the disturbance of the applied electric current density by an eigenflux field generated by inclusion with a proper transformation electric field.

Based on this method, the disturbed current density inside the matrix and particle domains are expressed as:

$$J_i-J_i^0=K^m{E}_i^d\mathrm{in}{\mathfrak{D}}^m,$$

(9a)

$$J_i-J_i^0=K^m[E_i^d-E_i^{*}]\mathrm{in}{\mathfrak{D}}^p.$$

(9b)

Here, we introduce $E_i^{*}$ , the transformation electric field due to the uniformly distributed electric field. The second-order tensor $S_{\mathrm{ij}}$ is analogous to the so-called Eshelby tensor in elasticity and depends only on the geometry of the ellipsoidal inclusion. It connects $E_i^d$ , and $E_i^{*}$ inside the inclusion domain as:

$$E_i^d=S_{\mathrm{ij}}{E}_j^{*}\mathrm{in}{\mathfrak{D}}^p.$$

(10)

The formulas to calculate $S_{\mathrm{ij}}$ for inclusions with different geometries can be found in Appendix 1. The resultant electric current density in the particle is given by the sum of equations (8) and (9b) as:

$$J_i=K^m[E_i^0+E_i^d-E_i^{*}]\mathrm{in}{\mathfrak{D}}^p.$$

(11)

Upon comparing equations (7) and (11), we get:

$$K^m[E_i^0+E_i^d-E_i^{*}]=K_{\mathrm{ij}}^p[E_j^0+E_j^d]\mathrm{in}{\mathfrak{D}}^p.$$

(12)

The above equation expresses a relation between the real composite and the equivalent inclusion.

2.3. Matrix with numerous ellipsoidal particles

The approach presented in section 2.2 for a single ellipsoidal particle can be extended for an infinite matrix with numerous ellipsoidal particles as shown in Figure 2. In this case, the interactions between the particles should also be taken into account.

Figure 2.

Numerous ellipsoidal particles embedded in an infinite matrix.

Here, the total electric field $E_i$ is decomposed into three parts as:

$$E_i=E_i^0+E_i^d+E_i^n\mathrm{in}\mathfrak{D},$$

(13)

where $E_i^n$ is the average disturbance of the electric field in the matrix due to conductive networks of particles and is defined as:

$$E_i^n=\frac{1}{V_{{\mathfrak{D}}^m}}{\int }_{{\mathfrak{D}}^m}[E_i-E_i^0]\mathrm{d}V.$$

(14)

With a similar approach to the EIM for a single ellipsoid presented in the previous section, equations (9b) and (12) can be rewritten for composites with numerous ellipsoidal particles as:

$$J_i-J_i^0=K^m[E_i^d+E_i^n-E_i^{*}]\mathrm{in}{\mathfrak{D}}^p,$$

(15)

and

$$K^m[E_i^0+E_i^d+E_i^n-E_i^{*}]=K_{\mathrm{ij}}^p[E_j^0+E_j^d+E_j^n]\mathrm{in}{\mathfrak{D}}^p.$$

(16)

By defining |□| that specify the outer boundary of □ and the partial differentiation operator ${\partial }_i\square =\partial \square /\partial x_i$ , and integrating over the entire domain $\mathfrak{D}$ , the electric current density $J_i-J_i^0$ vanishes using Gauss’ divergence theorem as:

$${\int }_{\mathfrak{D}}[J_i-J_i^0]\mathrm{d}V={\int }_{|\mathfrak{D}|}[J_j-J_j^0]n_j{x}_i\mathrm{d}V-{\int }_{\mathfrak{D}}{\partial }_j[J_j-J_j^0]x_i\mathrm{d}V=0,$$

(17)

since $(J_j-J_j^0)n_j=0$ on $|\mathfrak{D}|$ and ${\partial }_j(J_j-J_j^0)=0$ in $\mathfrak{D}$ . Hence, by integrating over the entire domain of the composite, equation (15) becomes:

$${\int }_{\mathfrak{D}}{E}_i^n\mathrm{d}V={\int }_{\mathfrak{D}}[E_i^{*}-E_i^d]\mathrm{d}V,$$

(18)

that results in:

$$E_i^n=\frac{1}{V_{\mathfrak{D}}}{\int }_{{\mathfrak{D}}^p}[E_i^{*}-E_i^d]\mathrm{d}V.$$

(19)

Thus, recalling equation (13), the volumetric average of the total electric field is given as:

$$\overline{E_i}=\frac{1}{V_{\mathfrak{D}}}{\int }_{\mathfrak{D}}{E}_i\mathrm{d}V=E_i^0+\frac{1}{V_{\mathfrak{D}}}{\int }_{\mathfrak{D}}[E_i^n+E_i^d]\mathrm{d}V.$$

(20)

Making use of equations (15) and (17), we get:

$$\overline{E_i}=E_i^0+\frac{1}{V_{\mathfrak{D}}}{\int }_{{\mathfrak{D}}^p}{E}_i^{*}\mathrm{d}V.$$

(21)

The electric field when integrated over the inclusions can be written as:

$${\int }_{{\mathfrak{D}}^p}{E}_i\mathrm{d}V=K^m{[K^m{\delta }_{\mathrm{ij}}-K_{\mathrm{ij}}^p]}^{-1}{\int }_{{\mathfrak{D}}^p}{E}_j^{*}\mathrm{d}V.$$

(22)

Hence, using $\overline{E_i}$ from equation (21) and ${\int }_{{\mathfrak{D}}^p}{E}_i\mathrm{d}V$ from equation (22), the non-zero components of overall electrical conductivity $K_{\mathrm{ij}}$ in equation (6) can be estimated as:

$$K_{11}=K_{22}=K^m[1-\frac{{\int }_{{\mathfrak{D}}^p}{E}_1^{*}\mathrm{d}V}{V_{\mathfrak{D}}{E}_1^0+{\int }_{{\mathfrak{D}}^p}{E}_1^{*}\mathrm{d}V}],$$

(23a)

$$K_{33}=K^m[1-\frac{{\int }_{{\mathfrak{D}}^p}{E}_3^{*}\mathrm{d}V}{V_{\mathfrak{D}}{E}_3^0+{\int }_{{\mathfrak{D}}^p}{E}_3^{*}\mathrm{d}V}].$$

(23b)

Next, we need to determine ${\int }_{{\mathfrak{D}}^p}{E}_i^{*}\mathrm{d}V$ in the above expression based on the distribution of particles.

2.4. Composites with randomly distributed inclusions

In order to evaluate the overall electrical conductivity of a composite reinforced by randomly distributed inclusions, we need to determine ${\int }_{{\mathfrak{D}}^p}{E}_i^{*}\mathrm{d}V$ . The model used by Takao et al. [53] for the composite reinforced by randomly distributed inclusions is used in this paper. The orientation of every inclusion is described by Euler angles $\theta$ and $\psi$ as shown in Figure 2. Using the local coordinates $\{{\hat{x}}_1,{\hat{x}}_2,{\hat{x}}_3,\}$ , where the ${\hat{x}}_3$ -axis is set to coincide with the inclusion axis, the EIM from equation (16) for a representative inclusion yields:

$$K^m[{\hat{E}}_i^0+{\hat{E}}_i^d+{\hat{E}}_i^n-{\hat{E}}_i^{*}]=K_{\mathrm{ij}}^p[{\hat{E}}_j^0+{\hat{E}}_j^d+{\hat{E}}_j^n]\mathrm{in}{\mathfrak{D}}^p,$$

(24)

and

$${\hat{E}}_i^d=S_{\mathrm{ij}}{\hat{E}}_j^{*}\mathrm{in}{\mathfrak{D}}^p,$$

(25)

where the hat accent refers to the local coordinate system on representative inclusion. Equations (24) and (25) can be held for every inclusion distributed in the matrix. From equations (24) and (25), we can obtain:

$$[K^m{\delta }_{\mathrm{ij}}-K_{\mathrm{ij}}^p][{\hat{E}}_j^0+{\hat{E}}_j^n]=[K_{\mathrm{ij}}^p-K^m{\delta }_{\mathrm{ij}}]S_{\mathrm{jk}}{\hat{E}}_k^{*}+K^m{\hat{E}}_i^{*}\mathrm{in}{\mathfrak{D}}^p.$$

(26)

Upon introduction of second-order tensors $A_{\mathrm{ij}}$ and $B_{\mathrm{ij}}$ :

$$A_{\mathrm{ij}}=[K_{\mathrm{ik}}^p-K^m{\delta }_{\mathrm{ik}}]S_{\mathrm{kj}}+K^m{\delta }_{\mathrm{ij}},$$

(27a)

$$B_{\mathrm{ij}}=A_{\mathrm{ik}}^{-1}[K^m{\delta }_{\mathrm{kj}}-K_{\mathrm{kj}}^p],$$

(27b)

the value of ${\hat{E}}_i^{*}$ in equation (26) can be stated in a more compact way as:

$${\hat{E}}_i^{*}=B_{\mathrm{ij}}[{\hat{E}}_j^0+{\hat{E}}_j^n].$$

(28)

From equations (25) and (28), we can eliminate ${\hat{E}}_i^{*}$ and obtain:

$${\hat{E}}_i^d=S_{\mathrm{ij}}{B}_{\mathrm{jk}}[{\hat{E}}_k^0+{\hat{E}}_k^n].$$

(29)

The rotation tensor $Q_{\mathrm{ij}}$ defined as $a_i=Q_{\mathrm{ij}}{\hat{a}}_j$ links the local coordinate system of the particle and the global coordinate system in the matrix. This second-order tensor can be written in a matrix form as:

$$[Q_{\mathrm{ij}}]=\left[\begin{array}{ccc}\cos \psi \cos \theta \hfill & -\sin \psi \hfill & \cos \psi \sin \theta \hfill \\ \multicolumn{1}{c}{\cos \theta \sin \psi }& \cos \psi \hfill & \sin \psi \sin \theta \hfill \\ \multicolumn{1}{c}{-\sin \theta }& 0\hfill & \cos \theta \hfill \end{array}\right].$$

(30)

Using this tensor, we can transform the conductivity tensor into the global coordinate system. So, by transforming equations (28) and (29), we can obtain:

$$E_i^{*}=Q_{\mathrm{ij}}{B}_{\mathrm{jk}}{Q}_{\mathrm{kl}}^{-1}[E_l^0+E_l^n],$$

(31a)

$$E_i^d=Q_{\mathrm{ij}}{S}_{\mathrm{jk}}{B}_{\mathrm{kl}}{Q}_{\mathrm{lo}}^{-1}[E_o^0+E_o^n].$$

(31b)

Now, consider a unit sphere where the orientation of inclusion is indicated by Euler angles $\theta$ and $\psi$ . The integration of quantity of $E_i^{*}$ over particle domain ${\mathfrak{D}}^p$ is stated as:

$${\int }_{{\mathfrak{D}}^p}{E}_i^{*}\mathrm{d}V=\int \int E_i^{*}\xi (\theta ,\psi )\mathrm{d}s=\underset{\psi =0}{\overset{2\pi }{\int }}\underset{\theta =0}{\overset{\alpha }{\int }}{E}_i^{*}\xi (\theta ,\psi )\sin \theta \mathrm{d}\theta \mathrm{d}\psi ,$$

(32)

where $\alpha$ is the limit of the inclusion orientation angle $\theta$ as $|\theta |<\alpha$ and $\alpha \in [0,\pi /2]$ . The orientation distribution function $\xi (\theta ,\psi )$ is defined as the number of particles $(\rho )$ intersecting a unit area of the unit sphere multiplied by the volume of a single particle $V_{\mathrm{sp}}$ , i.e., $\xi (\theta ,\psi )=\rho V_{\mathrm{sp}}$ . The function $\xi$ is identity for a uniform random distribution but can take more general forms as shown in further sections. Similarly, we can compute the integration of the quantity of disturbed electric field over the particle domain as:

$${\int }_{{\mathfrak{D}}^p}{E}_i^d\mathrm{d}V={\int }_0^{2\pi }{\int }_0^{\alpha }{E}_i^d\xi (\theta ,\psi )\sin \theta \mathrm{d}\theta \mathrm{d}\psi .$$

(33)

The volume fraction of inclusions, $\varphi$ , is evaluated as:

$$\varphi =\frac{1}{V_{\mathfrak{D}}}{\int }_{{\mathfrak{D}}^p}\xi (\theta ,\psi )\mathrm{d}s=\frac{1}{V_{\mathfrak{D}}}{\int }_0^{2\pi }{\int }_0^{\alpha }\xi (\theta ,\psi )\sin \theta \mathrm{d}\theta \mathrm{d}\psi .$$

(34)

By inserting $E_i^{*}$ in equation (31a) into equation (32), and inserting the expression for $E_i^d$ from equation (31b) into equation (33), we arrive at:

$${\int }_{{\mathfrak{D}}^p}{E}_i^{*}\mathrm{d}V=\left[{\int }_0^{2\pi }{\int }_0^{\alpha }\xi \left[Q_{\mathrm{ik}}{B}_{\mathrm{kl}}{Q}_{\mathrm{lj}}^{-1}\right]\sin \theta \mathrm{d}\theta \mathrm{d}\psi \right][E_j^0+E_j^n],$$

(35a)

$${\int }_{{\mathfrak{D}}^p}{E}_i^d\mathrm{d}V=\left[{\int }_0^{2\pi }{\int }_0^{\alpha }\xi \left[Q_{\mathrm{ik}}{S}_{\mathrm{kl}}{B}_{\mathrm{lm}}{Q}_{\mathrm{mj}}^{-1}\right]\sin \theta \mathrm{d}\theta \mathrm{d}\psi \right][E_j^0+E_j^n].$$

(35b)

Substituting $V_{\mathfrak{D}}$ from equation (34) into equation (19) and introducing tensors $C_{\mathrm{ij}}$ and $D_{\mathrm{ij}}$ as:

$$C_{\mathrm{ij}}=\frac{{\int }_0^{2\pi }{\int }_0^{\alpha }\xi \left[Q_{\mathrm{ik}}{B}_{\mathrm{kl}}{Q}_{\mathrm{lj}}^{-1}\right]\sin \theta \mathrm{d}\theta \mathrm{d}\psi }{{\int }_0^{2\pi }{\int }_0^{\alpha }\xi \sin \theta \mathrm{d}\theta \mathrm{d}\psi },$$

(36a)

$$D_{\mathrm{ij}}=\frac{{\int }_0^{2\pi }{\int }_0^{\alpha }\xi \left[Q_{\mathrm{ik}}{S}_{\mathrm{kl}}{B}_{\mathrm{lm}}{Q}_{\mathrm{mj}}^{-1}\right]\sin \theta \mathrm{d}\theta \mathrm{d}\psi }{{\int }_0^{2\pi }{\int }_0^{\alpha }\xi \sin \theta \mathrm{d}\theta \mathrm{d}\psi },$$

(36b)

to compact the formulation further, equation (19) yields:

$$E_i^n=\varphi [C_{\mathrm{ij}}-D_{\mathrm{ij}}][E_j^0+E_j^n].$$

(37)

Considering the fact that $C_{\mathrm{ij}}$ and $D_{\mathrm{ij}}$ are the diagonal matrices due to $K_{\mathrm{ij}}^p$ , $K^m{\delta }_{\mathrm{ij}}$ , and $S_{\mathrm{ij}}$ being diagonal, the components of $E_i^n$ can be explicitly expressed in terms of the components of vector $E_i^0$ as:

$$E_1^n=E_2^n=\frac{\varphi [C_{11}-D_{11}]E_1^0}{1-\varphi [C_{11}-D_{11}]},$$

(38a)

$$E_3^n=\frac{\varphi [C_{33}-D_{33}]E_3^0}{1-\varphi [C_{33}-D_{33}]},$$

(38b)

so that the components of ${\int }_{{\mathfrak{D}}^p}{E}_i^{*}\mathrm{d}V$ in equation (35a) are written as:

$${\int }_{{\mathfrak{D}}^p}{E}_1^{*}\mathrm{d}V={\int }_{{\mathfrak{D}}^p}{E}_2^{*}\mathrm{d}V=\frac{\varphi V_{\mathfrak{D}}{C}_{11}{E}_1^0}{1-\varphi [C_{11}-D_{11}]},$$

(39a)

$${\int }_{{\mathfrak{D}}^p}{E}_3^{*}\mathrm{d}V=\frac{\varphi V_{\mathfrak{D}}{C}_{33}{E}_3^0}{1-\varphi [C_{33}-D_{33}]}.$$

(39b)

Finally, the non-zero components of the overall electrical conductivity $K_{\mathrm{ij}}$ of the composite in equation (23) can be estimated as:

$$K_{11}=K_{22}=K^m[1-\frac{\varphi C_{11}}{1+\varphi D_{11}}],$$

(40a)

$$K_{33}=K^m[1-\frac{\varphi C_{33}}{1+\varphi D_{33}}].$$

(40b)

The established EIM for estimating the overall electrical conductivity of randomly distributed inclusions inside a matrix in equations (40) is only valid when there is no electron hopping between particles, and thus, no conductive networks of particles are formed in the composite.

2.5. Particles with interphase coatings

To take into account the effects of electron tunneling and conductive networks, the EIM derived in section 2.4 for randomly distributed particles is extended for a system with particles with interphase layers by replacing every regular particle with an equivalent particle. For this purpose, it is assumed that every particle is surrounded by an interphase layer of thickness $t$ and isotropic electrical conductivity of $K^{\mathrm{int}}$ as shown in Figure 3. Here, the particles are chosen to be cylinders.

Figure 3.

An equivalent cylinder particle with its interphase layer surrounding it divided into three parts.

The effective conductivity of the equivalent cylinder can be derived using Maxwell’s equations and the rule of mixture [41, 54]. This requires the equivalent cylinder to be divided into two isotropic interphases (parts I and III) and one transversely isotropic region of CNT/interphase (part II) as shown in Figure 3. The overall electrical conductivity of part II along ${\hat{x}}_3$ is denoted by $K_{33}^{\mathrm{II}}$ and can be evaluated as:

$$K_{33}^{\mathrm{II}}=\frac{K_{33}^p{D}^2+K^{\mathrm{int}}[4\mathrm{Dt}+4t^2]}{{[D+2t]}^2}.$$

(41)

The overall electrical conductivity of part II along ${\hat{x}}_1$ is denoted as $K_{11}^{\mathrm{II}}$ and is evaluated by applying a test electric field $E^t$ on the equivalent cylinder along the ${\hat{x}}_1$ -axis. Maxwell’s equations require the electric scalar potential $U$ to satisfy Poisson’s equation that is given in the cylindrical coordinate system $\{r,\theta ,z\}$ as:

$${\nabla }^{2}U=\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial U}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}U}{\partial {\theta }^2}=0.$$

(42)

The boundary conditions are prescribed as:

$$U^p{|}_{r=0}=\mathrm{constant},E^m{|}_{r\to \infty }=-\frac{\partial U^m}{\partial r}{|}_{r\to \infty }=E^t,$$

(43)

$$U^p{|}_{r=D/2}=U^{\mathrm{int}}{|}_{r=D/2},-K_{11}^p\frac{\partial U^p}{\partial r}{|}_{r=D/2}=-K^{\mathrm{int}}\frac{\partial U^{\mathrm{int}}}{\partial r}{|}_{r=D/2},$$

(44)

$$U^{\mathrm{int}}{|}_{r=D/2+t}=U^m{|}_{r=D/2+t},-K^{\mathrm{int}}\frac{\partial U^{\mathrm{int}}}{\partial r}{|}_{r=D/2+t}=-K^m\frac{\partial U^m}{\partial r}{|}_{r=D/2+t},$$

(45)

where the superscripts $p,m$ , and $\mathrm{int}$ refer to particle, matrix, and interphase, respectively. The above set of equations for $U$ results in:

$$U^p=2\gamma K^{\mathrm{int}}{E}^{t}r\cos \theta ,0\le r\le D/2,$$

(46)

$$U^{\mathrm{int}}=\gamma [K^{\mathrm{int}}+K_{11}^p+{\left[\frac{D}{2r}\right]}^2[K^{\mathrm{int}}-K_{11}^p]]E^{t}r\cos \theta ,D/2<r<D/2+t,$$

(47)

$$U^m=[\gamma [{\left[\frac{D+2t}{2r}\right]}^2[K^{\mathrm{int}}+K_{11}^p]+{\left[\frac{D}{2r}\right]}^2[K^{\mathrm{int}}-K_{11}^p]]+{\left[\frac{D+2t}{2r}\right]}^2-1]E^{t}r\cos \theta ,D/2+t<\infty ,$$

(48)

where:

$$\gamma =2K^m{[[K^{\mathrm{int}}-K_{11}^p][K^{\mathrm{int}}-K^m]{\left[\frac{D}{D+2t}\right]}^2-[K^{\mathrm{int}}+K_{11}^p][K^{\mathrm{int}}+K^m]]}^{-1}.$$

(49)

Also, the transverse electric field of the interphase and CNT parts of the equivalent cylinder can be evaluated as:

$$E_1^p=-\frac{\partial U^p}{\partial {\hat{x}}_1}=-\frac{\partial U^p}{\partial r}\frac{\partial r}{\partial {\hat{x}}_1}=-\frac{1}{\cos \theta }\frac{\partial U^p}{\partial r},$$

(50)

$$E_1^{\mathrm{int}}=-\frac{\partial U^{\mathrm{int}}}{\partial {\hat{x}}_1}=-\frac{\partial U^{\mathrm{int}}}{\partial r}\frac{\partial r}{\partial {\hat{x}}_1}=-\frac{1}{\cos \theta }\frac{\partial U^{\mathrm{int}}}{\partial r},$$

(51)

with:

$$J_1^p=K_{11}^p{E}_1^p,\mathrm{and}{J}_1^{\mathrm{int}}=K^{\mathrm{int}}{E}_1^{\mathrm{int}}.$$

(52)

Using the volumetric average, we can obtain:

$$\frac{1}{V}{\int }_V{J}_1\mathrm{d}V=K_{11}^{\mathrm{II}}\frac{1}{V}{\int }_V{E}_1\mathrm{d}V.$$

(53)

Therefore, using equations (52) and (53), we can write:

$$K_{11}^{\mathrm{II}}=\frac{D^2{K}_{11}^p{K}^{\mathrm{int}}+2K^{\mathrm{int}}[K^{\mathrm{int}}+K_{11}^p][t^2+\mathrm{Dt}]}{D^2{K}^{\mathrm{int}}+2[K^{\mathrm{int}}+K_{11}^p][t^2+\mathrm{Dt}]}.$$

(54)

Finally, the effective conductivity of an equivalent cylinder with a length of $l$ and diameter of $D$ along ${\hat{x}}_1,{\hat{x}}_2$ , and ${\hat{x}}_3$ is derived by:

$${\tilde{K}}_{11}^p={\tilde{K}}_{22}^p=\frac{K^{\mathrm{int}}}{l+2t}[2t+\frac{l{K}_{11}^p{D}^2/2+l[K_{11}^p+K^{\mathrm{int}}][t^2+\mathrm{Dt}]}{D^2{K}^{\mathrm{int}}/2+[K_{11}^p+K^{\mathrm{int}}][t^2+\mathrm{Dt}]}],$$

(55a)

$${\tilde{K}}_{33}^p=\frac{K^{\mathrm{int}}[l+2t][K_{33}^p{D}^2/4+K^{\mathrm{int}}[\mathrm{Dt}+t^2]]}{t{K}_{33}^p{D}^2/2+2t{K}^{\mathrm{int}}[\mathrm{Dt}+t^2]+l{K}^{\mathrm{int}}{[D/2+t]}^2},$$

(55b)

where the electrical conductivity tensor in the local coordinate system of the equivalent cylinder is given as $\mathrm{diag}({\tilde{K}}_{11}^p,{\tilde{K}}_{11}^p,{\tilde{K}}_{33}^p)$ .

Based on these effective cylinders, the effective volume fraction $\tilde{\varphi }$ can be calculated from the ratio of the volume of one equivalent particle $V^{\mathrm{ep}}$ to the volume of one particle $V^p$ multiplied by the volume fraction of original particles as:

$$\tilde{\varphi }=\frac{V^{\mathrm{ep}}}{V^p}\varphi =\frac{\pi {[D/2+t]}^2[l+2t]}{\pi {[D/2]}^{2}l}\varphi =\frac{4\varphi {[D/2+t]}^2[l+2t]}{l{D}^2}.$$

(56)

Therefore, the EIM described in section 2.4 to estimate the overall electrical conductivity of randomly distributed inclusion composites can be extended for composites with interphase layer by replacing the inclusions by equivalent particles, i.e., replacing $\varphi$ in equation (40) by $\tilde{\varphi }$ in equation (56) and replacing $K_{\mathrm{ij}}^p$ in equation (27) by ${\tilde{K}}_{\mathrm{ij}}^p$ in equation (55).

Remark. The established EIM in section 2 for estimating the overall electrical conductivity of randomly distributed inclusions inside a matrix with or without the interphase layer is equivalent to EIM for estimating the thermal conductivity of composite with or without the interphase layer. In this case, the electrical conductivity $K_{\mathrm{ij}}$ is replaced by the thermal conductivity ${\kappa }_{\mathrm{ij}}$ , the electric field $E_i$ is replaced by the gradient of temperature ${\partial }_{i}T$ , and the electric current density $J_i$ is replaced by the heat flux $q_i$ .

3. Electron tunneling and conductive networks

From this section onward, we focus on CNT/polymer nanocomposites and the inclusions are chosen to be CNTs. When the distance between two CNTs is less than a threshold value, electrons can transport from one CNT to another via a phenomenon commonly referred to as electron tunneling or electron hopping [55–58]. In this section, we specialize the EIM developed for composites with interphase layers around inclusions in section 2 by considering the interphase layer around CNTs to be a model for electron hopping [38].

Simmons [59] derived a generalized method to calculate the tunneling current between electrodes separated by an insulating region. By assuming a uniform thickness of the insulator and neglecting any variations in barrier height along its thickness, we can apply the formula designed for a rectangular potential barrier. Based on this, a generalized framework is developed where the fillers are treated as three-dimensional (3D) continuum objects. This approach is widely used in the literature to model the conductivity due to electron tunneling for CNT-based composites [12, 41, 43, 44, 49, 60, 61]. The tunneling resistance $R_t$ at a junction between two CNTs is estimated as:

$$R_t=\frac{h^{2}d}{A_t{e}^2\sqrt{2m\Delta E}}\exp \left(\frac{4\pi d}{h}\sqrt{2m\Delta E}\right),$$

(57)

where $h$ is Planck’s constant, $d$ is the distance between two adjacent CNTs, $A_t$ is the area available for tunneling, $e$ is the charge of an electron, $m$ is the mass of an electron, and $\Delta E$ is the energy barrier which is equal to the work required for an electron to tunnel.

The distance $d$ between CNTs is chosen to be different for the electron hopping mechanism in comparison to the mechanism based on conductive networks. In order to capture the effect of electron tunneling, the thickness of the interphase layer is considered as $t=d_c/2$ , while the conductivity of the layer is described as $K^{\mathrm{int}}=d_c/(A_t{R}_t)$ . In the next section, the value of the gap distance between CNTs $d$ , is determined for conductive networks of CNTs.

3.1. Conductive networks

In order to consider the effect of the conductive networks of CNTs, the average junction distance between CNTs $d_a$ is introduced so that the thickness of the interphase layer is assumed as $t=d_a/2$ , while the conductivity of the layer is given by $K^{\mathrm{int}}=d_a/(A_t{R}_t)$ . The distance $d_a$ is estimated by a power-law as [41, 62]:

$$d_a={\left[\frac{{\varphi }_c}{\varphi }\right]}^{1/3}{d}_c,$$

(58)

where ${\varphi }_c$ is the critical volume fraction of percolation (percolation threshold), which is the starting volume fraction for the formation of the continuous conductive networks.

Eventually, the overall electrical conductivity of CNT-reinforced polymer nanocomposite can be estimated based on a simple rule of mixture considering the effect of electron hopping and conductive networks of CNTs. Estimating the fraction of CNTs taking part in forming conductive networks after percolation threshold as [39]:

$$\zeta =\frac{{\varphi }^{1/3}-{\varphi }_c^{1/3}}{1-{\varphi }_c^{1/3}};({\varphi }_c\le \varphi <1),$$

(59)

while the rest contribute to the electron tunneling mechanism, the overall electrical conductivity of the composite is:

$$K=(1-\zeta )K_{\mathrm{T}}+\zeta K_{\mathrm{N}}.$$

(60)

here, $K_T$ is the contribution to the electrical conductivity from the tunneling effect which is evaluated by the extended EIM derived in section 2 by substituting the thickness $t=d_c/2$ and conductivity $K^{\mathrm{int}}=d_c/(A_t{R}_t)$ of interphase in equation (55). $K_N$ is the contribution to the electrical conductivity from conductive networks and is evaluated by the extended EIM derived in section 2 by substituting the thickness $t=d_a/2$ and conductivity $K^{\mathrm{int}}=d_a/(A_t{R}_t)$ of interphase in equation (55) under the assumption $L\to \infty$ in equation (10).

It is noted that we started from a general equivalent inclusion model with a wider range of applications in section 2 which is valid in predicting the thermal or electrical conductivity of any type of composite reinforced by particles that can be assumed as ellipsoids with various aspect ratios. In section 2.5, the general EIM was specified for fibers (cylinder-shaped fillers) to consider the effect of the interphase layer. It is noteworthy that this interphase layer consideration for fibers can be re-evaluated when dealing with other shapes of fillers. The model did not account for quantum phenomena such as electron hopping and conductive networks of particles yet. The effects of electron tunneling and conductive networks were considered into account in section 3 for CNT/polymer composites. It should be noted that this electron tunneling and conductive networks consideration for CNTs can be re-evaluated when dealing with other nanoparticles. Analytical models have some assumptions which simplify the actual problem in the real world. It is acknowledged that the presented model comes with certain limitations. For example, it does not account for commonly encountered features such as CNT waviness, agglomeration, segregation, variability in dimensions, multiple coatings, and possibilities of void formation [12, 43, 44, 60].

4. Results and discussion

In this section, the model developed in section 3 is validated using experimental data and is used to study the influence of CNT dimension, volume fraction, orientation, and distribution on the overall conductivity of the nanocomposite. The aspect ratio of CNTs is considered to be high $(l/D>>1)$ , and all CNTs are assumed to be isotropic $(K_{11}^p=K_{33}^p)$ [12, 38]. A value of $d_c=1.8$ nm is used for the cut-off distance between CNTs to allow for electron-hopping as has been reported in several papers [41, 63, 64]. The parameters $h$ , $e$ , and $m$ are the universal constants. The values of these parameters are listed in Table 1.

Table 1.

Model parameters with constant values.

Parameter	Value	Units
$d_c$	1.8	nm
h	$6.62607\times {10}^{-34}$	${\mathrm{m}}^2\mathrm{kg}/\mathrm{s}$
e	$-1.60218\times {10}^{-19}$	C
m	$9.10938\times {10}^{-31}$	kg

However, choosing other parameters needs more care since different types of CNTs and polymer matrices with various manufacturing processes for a wide range of applications have been reported. Accordingly, we estimate the model parameters by validation with several experimental studies in section 4.1. The length $l$ , diameter $D$ , and the percolation threshold ${\varphi }_c$ of CNTs are directly chosen from the measured or reported values in those works. Deng and Zheng [39] and Gao and Li [65], among others, have used an analytical expression to predict the percolation threshold:

$${\varphi }_c(H)=\frac{9H[1-H]}{2+15H-9H^2},$$

(61)

where $\gamma$ is the aspect ratio of the CNTs defined as $\gamma =l/D$ , and $H(\gamma )=\frac{1}{{\gamma }^2-1}[\frac{\gamma }{\sqrt{{\gamma }^2-1}}\ln (\gamma +\sqrt{{\gamma }^2-1})-1]$ . However, our investigation indicates that the above expression does not give an accurate estimation for every case study and as denoted in some papers in the literature, the value of ${\varphi }_c$ depends on the characteristics of polymer and CNTs and is hard to predict by an analytical expression for all types of polymer composites.

Nevertheless, the values of electrical conductivity of CNT nanoparticles $K_{11}^p$ and the energy barrier $\Delta E$ of various polymers are reported in a very wide range and the effective conductivity is very sensitive to these values. We determine these parameters separately for each experiment in section 4.1 by a least squares minimization technique. After model validation, parameters for the case of polystyrene (PS) reinforced by multi-walled CNTs [60] are used for the parametric study in section 4.2.

We utilize Mathematica software to perform numerical study cases in this study. The method proves highly efficient and economical in terms of time, thanks to the inherent advantages of this analytical approach. The Mathematica code for this analysis is shared in an open-source format [66]. The link to the code can be accessed in “Availability of data and materials” at the end of this manuscript.

4.1. Validation with experimental data

In order to verify the accuracy of the presented model, the analytical results are compared with some experimental studies herein. The case studies are chosen from four different nanocomposites—CNT/PS from Wang et al. [60], CNT/thermoplastic polyurethane (TPU) from Takeda et al. [43], CNT/epoxy from Kim et al. [67], and CNT/thermoplastic ultrahigh molecular weight polyethylene (UHMWPE) from Lisunova et al. [68]. The values of length and diameter of CNTs, percolation volume fraction, and intrinsic electrical conductivity of matrices for all four examples are listed in Table 2.

Table 2.

Parameters used for model validation in section 4.1.

Parameter	Units	PS/CNT	TPU/CNT	Epoxy/CNT	UHMWPE/CNT
l	$\mu \mathrm{m}$	3	30	30	10
D	nm	15	10	15	16
${\varphi }_c$	%	0.46	0.0723	0.02	0.05
$K^m$	S/m	${10}^{-11}$	${10}^{-11}$	${10}^{-11}$	${10}^{-13}$
Reference		Wang et al. [60]	Mora et al. [43]	Kim et al. [67]	Lisunova et al. [68]

PS: polystyrene; CNT: carbon nanotube; TPU: thermoplastic polyurethane; UHMWPE: ultrahigh molecular weight polyethylene.

The unspecified parameters are the intrinsic electrical conductivity of CNTs $K_{11}^p$ and the energy barrier $\Delta E$ . These parameters are determined based on the given ranges of $K_{11}^p=\{10-{10}^6\}\mathrm{S}/\mathrm{m}$ [41, 69, 70], and $\Delta E=\{0.1-5.0\}$ $e$ V [12, 62, 71], in comparison to the experimental results. Thus, the following function which expresses the least square logarithmic difference between the analytical value and the experimental data is introduced:

$$f(K_{11}^p,\Delta E)=\sum _{n=1}^m{[\underset{10}{\log }(K_{11}(K_{11}^p,\Delta E,{\varphi }_i))-\underset{10}{\log }(K_{11}^{\exp }({\varphi }_i))]}^2.$$

(62)

In equation (62), $m$ is the number of points in the experimental data for different volume fractions ${\varphi }_i$ and $K_{11}^{\exp }$ is the electrical conductivity measured in the experiments. The goal is to minimize this function to reduce the difference between the analytical model and experiments and thereby determine optimized values of $K_{11}^p$ and $\Delta E$ . Optimization using equation (62) results in a value of $K_{11}^p=1000\mathrm{S}/\mathrm{m}$ $\pm 10\%$ for all the four cases. Garcia-Macias et al. [12] reported that analytical results with electrical conductivity of CNT in the range of $\{100-1000\}\mathrm{S}/\mathrm{m}$ agree well with the experiments. To simplify the problem, we fix $K_{11}^p=1000\mathrm{S}/\mathrm{m}$ for all the cases, which is the most commonly used value in the literature and we only analyze the influence of the value of $\Delta E$ .

Our computation results in the value of $\Delta E$ as $0.7,1.2,2.7$ , and $1.1$ eV for CNT/PS, CNT/TPU, CNT/epoxy, and CNT/UHMWPE, respectively. In comparison, Wang et al. [60] used $K_{11}^p={10}^4\mathrm{S}/\mathrm{m}$ and $\Delta E=2.5\mathrm{eV}$ to model the experiments on CNT/PS polymer. These parameter values differ significantly from the optimized range and were not experimentally measured but taken from prior modeling data on CNT nanocomposites. On the contrary, Mora et al. [43] calculated $K_{11}^p$ as 987 S/m and used a prior published value of $\Delta E=1.5\mathrm{eV}$ for modeling the experiments on CNT/TPU. These parameter values are very close to the optimized range obtained herein.

Figure 4 illustrates the electrical conductivity of PS/CNT, TPU/CNT, epoxy/CNT, and UHMWPE/CNT nanocomposites predicted by the current model versus those given by experiments for different volume fractions of CNT. The graphs are plotted for 0.1, 2.5, and 5.0 eV, and the optimized value of $\Delta E$ for each case study.

Figure 4.

Comparing present model to experiments for (a) PS/CNT by Wang et al. [60], (b) TPU/CNT by Mora et al. [43], (c) epoxy/CNT by Kim et al. [67], and (d) UHMWPE/CNT by Lisunova et al. [68]. The black curves with the label $\Delta E(o)$ show the conductivity for optimized values of the energy barrier.

All graphs show a good agreement between the model and the experimental results. In Figure 4(a), it can be observed that for energy barriers higher or lower than $0.7$ , the results are far from the experimentally obtained values. The same trend can be observed for other cases which indicates that determining the energy barrier plays an important role in predicting the electrical conductivity by the present model. Also, as expected, increasing the value of the energy barrier drastically decreases the electrical conductivity after percolation. This effect is more pronounced for lower volume fractions of CNTs. Moreover, graphs with a higher energy barrier have a softer jump after percolation, while graphs with a lower energy barrier show a sharp change after percolation. It should be noted again that since the distribution of CNTs is uniform in all the cases considered here, the composites are isotropic, i.e., $K_{11}=K_{22}=K_{33}$ .

4.2. Effect of CNT orientation, distribution, and dimensions

Having validated our model with experimental data, we now study the effect of CNT alignment, distribution, and dimensions on the transverse and longitudinal effective electrical conductivity of the nanocomposites. For this parametric study, PS/CNT composite material from Wang et al. [60] is selected. The values of the parameters for PS/CNT composite are listed in Table 2. Also, the value of the energy barrier is obtained as $\Delta E=0.7\mathrm{eV}$ in the previous section. These values are used for all calculations in this section unless otherwise stated.

4.2.1. Limit angle

As discussed in section 1, non-uniform distribution of CNTs can be achieved by application of large strains on isotropic composites or directly manufacturing composites with a preferred CNT alignment. Fabricating devices with well-aligned CNTs is also desirable for better performance with lesser materials and smaller designs [72]. One of the parameters that control the alignment of CNT particles is the limit angle $\alpha$ introduced in section 2.4. It is reminded that $\alpha$ is the limit angle of $\theta$ , which is the angle made between a CNT and $x_3$ direction as it was shown in Figure 2. Hence, the probability of an inclusion with a tilt angle $\theta >\alpha$ is zero. Figure 5 shows the influence of volume percentage of CNT $\varphi$ on the transverse and longitudinal electrical conductivity of composites with different limit angles, $\alpha =\{\pi /32,\pi /16,\pi /8,\pi /4,\pi /2\}$ .

Figure 5.

The effect of different $\alpha$ angles on the (a) transverse and (b) longitudinal electrical conductivity of the composite.

As expected, upon reducing $\alpha$ while maintaining the volume fraction, the transverse electrical conductivity increases while the longitudinal electrical conductivity decreases. Since the longitudinal conductivity is driven by long CNTs while the transverse conductivity is driven by a combination of CNTs and matrix, the impact of limit angle on $K_{11}$ is significantly more pronounced than that on $K_{33}$ on a logarithmic scale. However, the increase of $K_{33}$ by a decrease of $\alpha$ is more evident here in the linear scale. Figure 6 depicts the effect of limit angle $\alpha$ on the transverse and longitudinal conductivities on a linear scale for $\varphi =\{1\%,1.5\%,2\%\}>{\varphi }_c$ volume fraction of CNTs. These figures provide an insight into the effects of varying the limit angle of inclusion orientation on transverse and longitudinal electrical conductivity.

Figure 6.

The effect of different $\alpha$ angles on the transverse (solid lines) and longitudinal (dashed lines) electrical conductivity of the composite.

Figure 6 confirms that decreasing the angle from ${90}^o$ which represents the full random distribution states to near $0^o$ , which is for the fully aligned CNTs along $x_3$ -axis, results in a drastic decrease in $K_{11}$ so that in the limit $\alpha \to 0$ , it approaches the matrix conductivity $K_{11}\to K^m$ . For all volume fractions, $K_{33}$ increases about three times as $\alpha$ drops from ${90}^o$ to $0^{°}$ . These results highlight that the effect of the limit angle on conductivity has a weak correlation with the volume fraction of CNT, particularly for small $\varphi$ . As anticipated, for a fully random distribution $(\alpha ={90}^o)$ , the transverse and longitudinal conductivities are equal $(K_{11}=K_{33})$ . This observation aligns with expectations, further affirming the validity of our model in capturing the conductivity behavior across various alignment scenarios.

4.2.2. Distribution function

The other parameter that directs the orientations of inclusions is the distribution function $\xi$ introduced in section 2.4. Unlike the angle $\alpha$ that acts as a sharp cut-off threshold and is hard to measure experimentally, the distribution function provides a smooth transition for the probability of the presence of CNTs along a particular direction and is, therefore, more amenable to modeling approaches. It is also challenging to accurately define a specific limit angle. Using a distribution function, on the contrary, is more feasible to account for the variability in inclusion alignment, thereby providing a more realistic and flexible framework for controlling the alignment of inclusions.

In the preceding analysis, a uniform distribution $(\xi =1)$ , was used. In this section, we study the effect of using different distribution functions including a linear function, $\xi =2-4\theta /\pi$ , a cosine function, $\xi =\pi /2\cos (\theta )$ , and an exponential function, $\xi =\pi /1.25\exp (-2{\theta }^2)$ . The probability distribution function is normalized as:

$${\int }_0^{\pi /2}\xi (\theta )\mathrm{d}\theta =\frac{\pi }{2},$$

(63)

and is demonstrated in Figure 7 for the uniform, linear, cosine, and exponential distributions. As shown, the exponential function gives the most aligned inclusions in one direction while in uniform distribution which is an ideal case, the inclusions are randomly distributed. We note that the probability distribution functions defined above are not unique. They can also be defined such that the function reaches zero smoothly. For example, the cosine distribution function $\xi =\pi /2\cos (\theta )$ can be replaced by $\xi =\cos (\theta \pi /2\alpha )$ for $|\theta |<\alpha$ .

Figure 7.

The uniform, linear, cosine, and exponential distribution functions versus $\theta$ .

Figure 8 shows the effect of the four different distribution functions on the transverse and longitudinal electrical conductivity of composite for different volume fractions. As expected, a notable trend is observed; transitioning from a uniform distribution function to a non-uniform function exerts a pronounced effect on electrical conductivity. Specifically, such a shift results in higher longitudinal conductivity and lower transverse conductivity. This disparity becomes particularly evident and significant as the volume content of CNTs is increased. This finding underscores the critical role played by the distribution function in shaping the electrical properties of CNT-based composites, highlighting its significance as a key parameter in the design and optimization of such materials.

Figure 8.

The effect of different distribution functions on the transverse (solid lines) and longitudinal (dashed lines) electrical conductivity of the composite.

4.2.3. CNT dimensions

Different single-walled and multi-walled CNTs with various architectures such as armchair, chiral, and zigzag can be found in various diameters and lengths, offering a versatile selection of composite fillers [73, 74]. Here, the effect of the size of CNT length and diameter on the electrical conductivity is investigated. The analysis of how varying the size of CNTs impacts the electrical conductivity contributes to a deeper understanding of the relationship between CNT dimensions and conductivity, guiding the selection of appropriate CNT fillers for various applications.

Figure 9 demonstrates the effect of the size of CNT particles with different volume fractions on the electrical conductivity of the composite with uniform distribution. It can be observed that longer CNTs demonstrate higher electrical conductivity compared to shorter ones. This can be attributed to the increased number of conductive pathways formed by longer CNTs, allowing for more electrons to pass through the material. Conversely, shorter CNTs exhibit lower conductivity due to the reduced number of conductive sites available. Also, keeping the volume fraction constant, composites made with CNTs of smaller diameters tend to exhibit higher effective electrical conductivity compared to composites made with CNTs of larger diameters. This is because of the higher surface area to volume ratio of CNTs with smaller diameters, which enables a greater number of conductive paths for electrons passing through the composite. Upon increasing the length or decreasing the diameter of the CNT beyond a certain threshold, no significant change in the conductivity is observed. At this point, the conductivity reaches a plateau as the dominant factors affecting conductivity, such as intrinsic material properties and the distribution of the CNTs, become more influential.

Figure 9.

The effect of (a) the length and (b) the diameter of CNTs on the electrical conductivity of composite with uniform distribution for $\varphi =\{1\%,1.5\%,2\%\}$ .

Figure 10 shows the impact of the size of CNT particles with different distribution functions on the electrical conductivity behavior of the composite material with $\varphi =1\%$ volume fraction of CNTs. The effect of different distribution functions on the electrical conductivity of the composite with different lengths and diameters of CNTs is notable. In particular, non-uniform distribution functions demonstrate a higher degree of sensitivity of longitudinal electrical conductivity to the changes in the dimension of the CNTs compared to the uniform function.

Figure 10.

The effect of (a) the length and (b) the diameter of CNTs on the transverse (solid lines) and longitudinal (dashed lines) electrical conductivity of composite with $\varphi =1\%$ volume fraction and different distribution functions.

Overall, these results highlight the complex interplay between CNT length, diameter, distribution, and volume fraction in determining the electrical conductivity of the composite material. This complexity highlights the need for a detailed understanding of the combined effects of these parameters. The findings emphasize the importance of carefully considering these parameters as they significantly impact the electrical conductivity and, consequently, the performance of such materials in various applications. The conductivity plays a crucial role in determining the overall performance of such materials across a spectrum of applications. It is evident that variations in CNT length, diameter, distribution, and volume fraction contribute to the diverse electrical properties observed in the composite. This insight emphasizes the necessity for a thoughtful approach to the design and engineering of CNT/polymer nanocomposites.

5. Conclusion

In this study, the overall electrical conductivity of composites with distributed inclusions is investigated with a specific focus on the impact of orientation and distribution of inclusion. An analytical model is presented based on Eshelby’s equivalent inclusion deriving mathematical relationships through a rigorous step-by-step process. Despite the popularity of EIM for predicting conductivity, a lack of detailed formulation derivation is sensed in recent literature. This paper addresses this gap and highlights the commonly overlooked non-uniform distribution of inclusions. The model is developed to predict the electrical conductivity of composites with random as well as preferred distribution of the fillers. This EIM formulation is valid for predicting the electrical and thermal conductivity of composites reinforced by randomly oriented ellipsoid fillers. The method is further extended by considering an interphase layer around the inclusion.

Subsequently, the method is extended for electrical conductivity of CNT-filled composites by assuming the interphase layers around the particles function as electron tunneling and conductive networks. The validity of the model is confirmed by comparing its results with experimental data from four different CNT-reinforced polymer matrices. The values of the energy barrier are determined for each case study based on a comparison to the experimental data, and it is observed that the energy barrier plays a crucial role in predicting electrical conductivity. It is worth noting that this issue has been rather avoided to be discussed in depth in the literature, and the sensitivity of analytical models to these parameters is neglected in most cases except for a few articles. Furthermore, a parametric study is conducted to investigate the effects of different aspects of the reinforcements on the transverse and longitudinal electrical properties. The effects of CNT orientation and distribution patterns on the electrical conductivity, highlighting how the alignment of CNTs influences overall properties, are investigated. Two measurements, namely, the limit angle of inclusion orientation and probability distribution function, are used to control the orientation. The limit angle of CNT orientation is varied, and it is found that decreasing the angle from a uniformly random distribution to a fully aligned state results in a drastic decrease in the transverse electrical conductivity. However, the longitudinal electrical conductivity shows less sensitivity to the angle variation. Also, it is revealed that distributing CNTs with non-uniform probability distribution functions has a noticeable impact on electrical conductivity. In particular, this shift leads to an increase in longitudinal conductivity and a decrease in transverse conductivity. This difference becomes more prominent when the volume content of CNTs is raised. In addition, the effects of particle dimensions are examined. It is observed that composites with shorter CNTs demonstrate lower electrical conductivity compared to composites with longer CNTs. This is due to the reduced number of conductive pathways available in shorter CNTs. CNTs with smaller diameters exhibit higher electrical conductivity of composite compared to those with larger diameters of CNTs. However, beyond a certain threshold, further modification of the length or the diameter of CNTs has a minimal impact on electrical conductivity as other factors become more influential.

The presented analytical model facilitates the design and optimization of composite materials for specific electrical conductivity requirements. This approach finds applications in fields such as flexible electronics, smart materials, energy storage, and aerospace, guiding the development of advanced materials. The presented mathematical model assumes simplifications such as straight CNT particles with uniform dimensions, neglecting their actual curviness, variability, agglomeration, segregation, and void formation effects observed in reality. This analytical framework provides a straightforward method for predicting the electrical conductivity of CNT-reinforced composites. The same framework can be used to model the thermal or electrical conductivity of any composite where inclusions might have multiple coatings. However, for modeling CNT nanocomposites more accurately, certain manufacturing-induced complexities need to be accounted for. Further research directions may include expanding the model to incorporate additional factors and parameters such as the waviness and the clustering of CNTs or performing the model on other fillers such as graphite, as well as incorporating experimental data from a wider range of composite materials. Modeling these important features is an important research direction that will be considered as an extension of this work.