MODELING AND QUANTITATIVE ANALYSIS OF CONNECTIVITY AND CONDUCTIVITY IN RANDOM NETWORKS OF NANOTUBES

4 Increased interest to the study of nanotubes and nanotube-based systems originates from the unique properties of these nanoscale objects. For example, electrical conductivity of the carbon nanotube at 300 K is about 106 S/m, which is just an order of magnitude lower than that of the silver one, but at the same time, thermal conductivity is nearly the same as thermal conductivity of diamond (~3000 W×K/m). One outstanding feature of CNTs that has to be mentioned is mechanical strength – from this viewpoint CNTs outperform even the best existing materials such as Kevlar. That said, even more valuable properties can be of use when one combines already exceptional characteristics of nanotubes with specific advantages of conventional materials.


Introduction
Nanotubes are tube-like structures typically having up to a few tens of nanometers in diameter and length of up to several hundreds of nanometers. The most prominent example of such structures is indeed carbon nanotubes, commonly designated as CNTs. Early observations of carbon nanotubes were reported in [1], although an intense research of CNTs took its start much later [2]. An individual CNT is basically a cylinder created by rolling a layer of graphite up. In case multiple graphite layers are rolled-up, they form a structure called a multiwall carbon nanotube (MWCNT). In this respect, CNTs are sometimes referred to as SWCNTs (single-walled carbon nanotubes).

MODELING AND QUANTITATIVE ANALYSIS OF CONNECTIVITY AND CONDUCTIVITY IN RANDOM NETWORKS OF NANOTUBES A . S t e l m a s h c h u k
Speaking of electrical transport in CNT-based systems networks is shown to be dominated by resistance at network junctions, which scale with the size of the interconnecting bundles.
Examples that come to mind are epoxy resins reinforced with carbon nanotubes. In particular, it has been demonstrated that the storage modulus of epoxy resin rises substantially with the addition of CNTs [3], thermal conductivity of such nanocomposites grows with the CNT concentration [4] and electrical resistivity of epoxy decreases from 10 14 Ω·m (neat epoxy) to 10 Ω·m after doping with 0.4 % of CNTs [5]. Another important class of materials is PE-DOT-polymers with the addition of SWCNTs or MWCNTs as a nanofiller. Such nanocomposites can serve a number of purposes in various scientific, engineering, technological and biomedical applications [6].
However, there is a number of experimental problems to overcome in order to take advantage of reinforcement effect: for instance, aggregation and entanglement of nanotubes even at the low level of loadings have to be addressed. On the other hand, since physical properties, in particular, electrical conductivity, of nanocomposites are dependent on many factors, significant efforts are required to study them experimentally. Therefore, a number of approaches were developed in order to simulate the processes in nanotube-based composites using computational experiments.
Numerical simulations are capable of providing important insights that are very hard to achieve experimentally. Still, careful planning and realization of a computer experiment are far from being simple procedures. The model for calculations should be reliable and at the same has to retain a reasonable level of complexity so that the necessary computational resources (time and processing power) can be accumulated.
In view of the above-mentioned, it is particularly important to construct a comprehensive model that will consider structure-property connections between nanotubes and their electrical conductivity as well as geometric characteristics, such as alignment and aspect ratio. This work presents an attempt to develop such computational approach by combining the three-dimensional network structure simulation with the electrical model of tunneling transport. A systematic investigation carried out within the suggested framework can be useful for designing novel various CNT-polymer composites for electronic applications.

Literature review and problem statement
Networks of carbon nanotubes have proven to be important for modern electronics [7]. Creating CNT-based devices for applications is a complex multistage process that requires not only sophisticated experimental procedures but also thoughtful planning [8]. The most appropriate way to design the nanocomposite device is to start with numerical investigations considering the random nature of networks formed by nanotubes. Then, basic features of the system have to be established and the influence of crucial model parameters on these features can be understood.
When a mixture of interconnected tubes is dealt with, the important phenomenon is percolation [9] that occurs in random systems [10] and depends on a variety of factors such as ordering of elements [11] or specific features of the host matrix [12]. The problem to be addressed here is how the conductivity of a nanotube network inside a composite behaves depending upon the concentration of nanotubes, their alignment and geometrical factor [13]. Modeling is also important in view of sensing applications [14].
There is a vast variety of nanostructured sensor designs in the form of one-, two-layered thick-film structures [15] with different geometry of multifunctional junction [16]. Physical parameters of such elements are stable in time [17] unlike conventional sensor ceramics, which is characterized by thermally-induced electronic relaxation [18]. Such degradation transformations in these materials are initially described by stretched or suppressed exponential relaxation functions [19]. In addition, there are specific effects (for example, aqueous-voltaic effect) in ceramic-based nanomaterials) [20] that need modeling and quantitative analysis. Influence of free volumes on the conductivity of functional glass and ceramics [21] is still insufficiently investigated, so simulations can be extremely helpful when such sensor engineering is considered.
Various modeling approaches have been tried to date. These include using analytical models to predict the macroscopic response of a nanocomposite [22] and application of neural networks to simulate the influence of nanotubes on the properties of the given material [23]. Computer experiments on the effect of CNT agglomeration on the conductivity of CNT-based polymer composite were performed in [24]. Some interesting ideas were developed in [25], where the authors proposed the polyhedral finite-elements model and simulate the mechanical properties of CNTs.
Ultimately, the complexity of emerging applications and, therefore, respective numerical models, requires powerful computing resources. That is why one should always look for optimization possibilities when proposing a numerical model, and high-performance computing systems are often used for the purpose of carrying out computer experiments on large systems [26].
Although, as illustrated above, extensive computational efforts were put into the study of CNT networks, still there are almost no systematic studies that use 3D statistical modeling and assume different processes of conduction junctions formed by adjacent tubes. To bridge this gap, we are proposing a framework for numerical modeling of nanotube networks in the insulating medium, which allows for parametrisation, visualisation of the results and use of different formalisms for the description of electrical conductivity mechanisms.

The aim and objectives of the study
The aim of this work is to examine the behavior of DC electrical conductivity in "insulator-CNTs" model composite systems.
To achieve this aim, the following objectives have been identified: 1. Development of a reliable yet simple computer model, which can be used to investigate the percolation threshold with respect to the parameters of the system under study.
2. Description of conductivity between individual nanotubes at the junction points followed by the calculation of the total electrical conductivity of the generated network using matrix algebra and a combination of Kirchhoff's law and Ohm's law. 3. Verification of the model, computer experiments on nanotube networks with different parameters and estimation of computational efficiency depending on the size of the simulated volume.

Numerical model of nanotube network
The nanotube network can be thought of as a set of conductive nanotubes (NT) randomly placed within a 3D volume box (Fig. 1). This box is very often called a representative volume element (RVE). Nanotubes are considered as objects with a hard body, so that they do not penetrate each other by volume. This approach is more complex to simulate and results in the increased resources for computer simulations. On the other hand, this approach represents the behavior of real-world nanotube networks better. Two types of electric contact between NTs are taken into account: direct contact and contact when NTs are put at a very close distance and undergo a tunneling effect, which exists between nanoscale objects. Two opposite sides of the RVE are considered to serve as electrodes, to which external electrical voltage is applied. In order to define the total conductivity of the composite system, the conductivity between these electrodes has to be defined. Sets of interconnected tubes form aggregations and the size of these aggregations is rapidly growing as the volume concentration (loading) of nanotubes increases. At some point, concentration is large enough for the conduction across two electrodes to appear. This critical point is called the percolation threshold. For large nanotube networks, the percolation threshold value approaches to some constant value.
According to [27,28], tunneling conductivity has a more significant influence on the resulting conductivity of the RVE, than the resistance of the nanotubes themselves. So, one should not neglect the effect of tunneling conductivity while performing computer simulations of the nanotube composite behavior.

Fig. 1. Nanotube network bordered with two electrodes
The process of simulating the composite system starts with a random filling of the RVE volume with NTs. NTs with the predefined length and diameter are generated and added to the system one by one until achieving the desired volume concentration of NTs in the RVE.
In the simulation model, the nanotube is defined as a cylinder with spherical endings with the length l NT and diameter d. The central axis of the NT starts at the point A with the coordinates (x 1 , y 1 , z 1 ) and ends at the point B (x 2 , y 2 , z 2 ). The process of generating and random placing of a nanotube in the RVE consists of several phases. Randomly selected inside the RVE point A with the coordinates (x 1 , y 1 , z 1 ) stands for the start point of the cylinder axis [29]: where rand is a random floating point value within [0,1] range. Uniform distribution of NTs inside the RVE is expected to be achieved in this case. Then the random direction in space is described by two angles α and β: Based on these two direction angles, the position of the end point B is defined as: If needed, we cut off the exceeding part of an NT, so that we can be sure that every NT is actually located inside the RVE.
When generating a new NT and adding it into the system, we have to perform an additional check to verify that the newly generated tube does not overlap the existing ones. To perform this check, the shortest distance between the generated tube and all other NTs is calculated. Only in the case when all these distances are larger than the CNT diameter d, the newly generated NT is going to be added to the system. The next step after finishing the NTs generation is to find connections between them. Electrical contact between two NTs is formed if the shortest distance between them is shorter than the tunnel cut-off distance. Contact place between two NTs is defined by a pair of "junction" points: one on the axis of the first NT and the other one on the axis of the second NT. The process which yields both the shortest distance between NTs and the actual coordinates of these "junction" points is described below.
The coordinates of the arbitrary point on the NT axis can be described as follows: where p is the starting point of the NT, d  is the directional vector of the NT and α is a variable coefficient (α∈[0,1]), whose exact values define the points located on the NT axis.
Let us assume that p 1 +α· 1 d  and p 2 +β· 2 d  are the arbitrary points on the two NTs, respectively, and one needs to find the distance between them. Then (p 1 is the vector, which connects these two points. The minimal distance between NTs can be obtained as the norm of this aforementioned vector: We have to find the values of the coefficients α and β such that D has the minimal possible value. So, after a number of transformations, we obtain the following expressions for finding the α and β coefficients.
Auxiliary notations: The values of the unknown coefficients: If the obtained values of the coefficients α and β do not fall within the interval [0, 1], then their values should be adjusted to the nearest admissible values. The actual coordinates of "junction" points are calculated from the equation (9).
In order to accurately represent real NTs with the hardcore body, we need to introduce some additional restriction for the generated tubes. These restrictions state that NTs cannot overlap with each other by volume and the minimal distance between their surfaces is larger than the van der Waals separation distance. That is why we check the shortest distance between the central axis of the newly generated NT and the axis of all the other nanotubes already placed inside the RVE. If this shortest distance is smaller than the nanotube diameter plus the van der Waals separation distance, then it means that collision between those tubes exists and the newly generated tube has to be thrown away from the system. After all tubes are placed inside the RVE, the next problem is to locate interconnected tubes. We assume that the tubes are connected, if the closest distance between their surfaces is smaller than some tunneling cut-off distance. According to experimental results of measuring the tunneling effect conductivity between multiwalled NTs [30], this cutoff distance is around 2 nm.
In order to find out which one of the NTs actually takes part in the process of transmitting the current between two opposite sides of the nanocomposite sample, we utilized percolation theory approaches. A bunch of mutually connected nanotubes can form a cluster. The cluster is called conductive (or percolative) if it contains nanotubes connected to both electrodes. The existence of conductive cluster in the RVE is extremely important, because it allows including the RVE into an electric circuit. It is worth noting that the same RVE can have a number of independent conductive clusters. Mathematically, the system is represented as a graph and the percolation search is the search of the connected component of the graph, which contains elements that touch the opposite edges of the box ("electrodes"). "Electrodes" are simulated as pseudo-tubes (having zero-size dimensions).
In order to find the conducting cluster, we have implemented the weighted quick union algorithm with path compression. Thus, we define all conductive clusters in the composite system as there can be several parallel conductive clusters, especially at concentrations slightly larger than the percolation threshold for this particular system. Tubes not included in any of such clusters have no influence on the total conductivity of the system, therefore we ignore them in further calculations.
We have utilized the random resistor network approach in order to calculate the total equivalent conductivity of the composite system. A set of interconnected NTs is represented as a set of resistors connected to a junction point. Resistors represent both the intrinsic conductivity of the nanotube segment and the resistance of the contact between tubes. The network's junctions are the equivalents of points of contact between nanotubes. The coordinates of such points are selected so that the distances between them are equal to the minimum distance between nanotubes.
The tunneling effect at the "junction" points causes a contact resistance between a pair of nanotubes (Fig. 2). Suppose, the shortest distance between a pair of nanotubes is d kp , where d kp is shorter than d cutoff . Then, the contact resistance can be evaluated using the Landauer-Büttiker formalism as [31][32][33][34][35] where h is the Planck's constant; P is the transmission probability for the electron to tunneling between NTs; N is the number of conduction channels, which is dimensionless and related to the diameter of an NT [36]; e is the charge of an electron; d vdw is the van der Waals separation distance [37,38], which limits the minimum distance between a pair of NTs; d tunnel is the tunneling characteristic length; m e is the mass of an electron; ΔE is the height of energy barrier [39].

Fig. 2. Example of a contact pattern between nanotubes
Let us select a pair of points on the axis of an NT and define a part of an NT restricted by these two points. Suppose, the length of this segment is l s . Then, the intrinsic resistance of the NT segment can be found according to the formula [28,29]: 2 4 , where σ CNT is the intrinsic electrical conductivity of the NT and d denotes the diameter of the NT.
The random resistor network is formed by finding all "contact" points between NTs and local conductivities between them. The RVE is placed in the electric circuit. Since our random resistor network is now a part of this electric circuit, it is crucial to calculate its total conductivity. A combination of the Kirchhoff's current law (KCL) and Ohm's law is used to accomplish this goal.
Consider some element E connecting the nodes i and j in a resistor network. According to the Ohm's law, the current which enters and leaves this resistor element satisfies the expression below: where V i and V j are electric potential values at the nodes i and j, I i and I j define an electric current flowing through the i and j points, correspondingly. The local conductivity matrix e ij K is determined as follows: According to the KCL, the following system of linear algebraic equations for the entire resistor network is formed: Here K represents the global conductivity matrix of the resistor network. V i is an unknown electrical potential at each node, I 1 and I N are the currents entering and leaving the system, respectively.
Element-wise conductivity matrix e ij K is utilized for constructing the total matrix K. Initially, the matrix K is a zero-valued matrix.
But in most practical usages, the values of I 1 and I N are unknown. Instead, often the voltage difference applied to the RVE is set up: For that reason, first and last equations in the system (19) can be omitted and a new system of linear equations of dimensions (N-2)×(N-2) can be formed.
Actually, the nanotube network model shows a specific feature: only a small amount of resistors connected to each node, because of the low typical concentration of NTs inside the RVE for the system under study. So, the resulting matrix K is sparse and symmetric. The system of equations (19) is stored in a compressed column format, and a specific solver shipped with the SuperLU library [40][41][42] is utilized and optimized for handling sparse matrices. After solving the system and obtaining the values of the V i vector, all currents flowing through the resistor elements are calculated as: where σ ij is the conductivity of the resistor element between the i and j nodes and V i , V j are electrical potentials in those nodes. All values of electric currents I i can be utilized for emphasizing a conductive path between electrodes and depicting the distribution of current flow through the NT network.
In order to calculate the total current flowing through the resistor network, the currents through all resistor elements which are directly connected to one of the electrodes are added up: where m is the number of resistors being direct neighbors of the electrode. Finally, the total equivalent electrical conductivity of the random resistor network is calculated using the Ohm's law as follows: The obtained value of σ total represents conductivity of the nanotube network sample.

Connectivity and conductivity of nanotube network simulation results
The main purpose of the first numerical experiment is to study connectivity among nanotubes and the required conditions for the formation of a conductive path in the 3D random nanotube network depending on the parameters of the system. This experiment is important for further investigation of conductivity of the network since it provides data for evaluating the correctness of the model and its software implementation.
The output of this numerical experiment is the value of the percolation threshold, defined as volume concentration of NTs inside the RVE at which the system becomes conductive. Geometrical parameters of NTs, such as length and diameter, can influence the percolation threshold value. In the first numerical experiment, the percolation threshold as a function of the aspect ratio l NT /d has been investigated. Dimensions of the RVE are set to L X =L y =L Z =1,000 nm, the diameter of the nanotube is set to 2nm and by changing the length of nanotubes, different aspect ratios are achieved. In order to diminish the influence of specific nanotube configuration in space on the calculated conductivity value, the simulations are repeated 100 times for different nanotube configurations for the same concentration and nanotube aspect ratio. The series of computer simulations are performed for different aspect ratio values ranging from 60 to 160. The calculated results are depicted in Fig. 3. In the second computer experiment, the electrical conductivity of the nanotube network is calculated. The goal of this experiment is to compare the obtained results with the experimental measurements of the conductivity of the nanocomposite consisted of carbon nanotubes and polymer dielectric matrix, which are available, for example, in [29]. Such comparison is performed because the nanotube network model can represent electrical properties of such nanocomposite. In order to simulate such composite, the geometrical and electric parameters of the simulated composite system are set up as follows.
The size of the RVE along the X, Y and Z axes are 1.2 µm×1 µm×100 nm, respectively. The diameter d of NTs is set to 2 nm, the length l NT is 200 nm, so the aspect ratio value l NT /d of the NT is 100. The intrinsic conductivity of the NT segment is 10 4 S/m. The tunnel cut off distance is set to 1.8 nm. The average value among these simulations is assumed as the final value of total conductivity of the composite system under study. The obtained results show that for the concentration of 0.01, which is above the percolation threshold for the composite with the aforementioned parameters, the conductivity is equal to 5.5 S/m. This obtained value coincides with the experimental results provided in [29].
The dependency of electrical conductivity of the nanotube network on the nanotube aspect ratio for the model with the aforementioned parameters is also investigated. The series of computer simulations are performed for different aspect ratio values ranging from 60 to 160. The results of computer experiments are depicted in Fig. 4.
Then the influence of the cut-off distance value on the characteristic of the dependency of the composite total conductivity on the nanotube volume fraction is studied. In order to perform such analysis, the conductivity dependency for systems with the tunnel cut-off distance set to 1.5 nm, 1.8 nm, 2.1 nm, 2.5 nm is calculated. The results of these simulations are illustrated in Fig 5. The length of the NT is also changed to 300 nm, so the aspect ratio value l/d of the NT is 150.
From Fig. 5, one can make a conclusion that the tunneling cut-off distance does not influence significantly the way how the system behaves depending on volume fraction changes. Additionally, the tunneling cut-off distance has an effect on the actual values of the electric conductivity of the system. The calculation time of conductivity of the nanotube network as a function of the system size is analyzed. Configuration of the system at which tests were performed was as follows: Intel Core i7-5500U (4 cores) CPU, 2.4 GHz, 8 Gb RAM; Ubuntu 16.10 operation system, Mono software platform.
Average calculation time results are shown in Fig. 6. As one can see from Fig. 6, the calculation time rapidly increases with adding more tubes to the network. The time complexity of simulations is nearly quadratic one. In order to improve the calculations performance, the parallel computations can be utilized. Such advanced computational techniques allow using our model for calculating the conductivity of bigger networks and improving the accuracy of nanocomposite conductivity evaluations.

Discussion of the results of nanotube network conductivity and connectivity simulation
The proposed model, methodology and respective numerical experiment allowed establishing the conditions for the formation of a conductive path in the three-dimensional network of nanotubes with random distribution.
The most important findings are: -dependence of the threshold in the nanotube network on the concentration of nanotubes in the modeled volume; -dependence of the percolation threshold in the nanotube network on the nanotube aspect ratio; -increase of the СNT-composite conductivity with the tunneling cut-off distance.
These findings would serve as a useful tool when designing nanocomposites, as the optimal concentration of CNTs and their aspect ratio can be estimated prior to fabrication of samples with tailored electrical properties. But efficient application of the developed approach to realistic systems would require a significant increase of the simulated box size (and, respectively, the number of elements in the volume). And the downside here is that, as confirmed by the computational efficiency estimation, the computational time exponentially increases with the size of the system under study (the number of tubes in the network). To overcome this, parallel computing algorithms should be applied.
The performed study is a logical continuation of the previous research efforts and the next step is to compare the results of the conducted simulation with the results of measurements on real CNT-polymer composites, e. g. CNT-reinforced epoxy resins. On the computational side, further improvement of the discussed model will require considering the specific alignment of nanotubes inside the simulated volume as this factor can have a significant influence on the percolation threshold of the system.

Conclusions
1. The framework for computational studies of nanotube-insulator model composites has been developed. Within this framework, simulations of electrical conductivity considering a large number of interconnected nanotubes across the boundary edges of the representative volume element can be performed.
2. Numerical description of conductivity processes in random nanotube networks has been done. It was demonstrated, that within the proposed model, the conductivity percolation threshold depends on the concentration of nanotubes, their aspect ratio and tunneling distance.
3. The algorithm of total conductivity calculations using element-wise conductivity matrices was implemented. The calculated value of electrical conductivity of the model CNT-insulator composite above the percolation threshold is 5.5 S/m, which is comparable to those, reported from various experiments (1-10 S/m).
4. Simulation of the aspect ratio change in the range from 60 to 160 indicated a substantial effect on the percolation threshold, lowering it from 0.9 at 60 to 0.5 at 180. On the other hand, variation of the tunneling cut-off distance has a rather moderate influence on the conductivity of the system.