Nanoscale Continuum Modelling of Carbon Nanotubes by Polyhedral Finite Elements

As the geometry of a cell of carbon nanotube is hexagonal, a new approach is presented in modelling of single-walled carbon nanotubes using polyhedral finite elements. Effect of varying length, diameter, and thickness of carbon nanotubes on Young’s modulus is studied. Both armchair and zigzag configurations are modelled and simulated in Mathematica. Results from current approach found good agreement with the other published data.


Introduction
Since carbon nanotubes (CNTs) were first discovered in 1991 [1], they have gained attention from researchers due to their extraordinary physical properties, two of which are high stiffness and low density, making them ideal to be utilized as fibers in nanocomposite materials.Many experimental works have been performed, in order to determine the physical properties of different types of CNTs.These experiments are carried out at the nanoscale, by utilizing transmission electron microscopy (TEM) and atomic force microscopy (AFM).Their results contain large variabilities, due to the technological difficulties associated with the nanoscale.Young's modulus for single-walled carbon nanotubes (SWCNTs) and multiwalled carbon nanotubes (MWCNTs) is reported to be within the vast range of 0.32 TPa to 4.15 TPa and 0.27 TPa to 460 TPa, respectively [2][3][4][5][6][7][8][9][10][11][12][13].
Later, computational methods were employed to determine properties of CNTs, due to the difficulties and high cost faced in the experimental analysis.Computational simulation of CNTs is categorized into three, which are atomistic modelling, continuum modelling, and nanoscale continuum modelling.Atomistic modelling involves prediction of motion of each atom due to the interatomic forces and surrounding boundary conditions.The collective behavior of the atoms helps to predict the material behavior, such as deformation, phase change, or other phenomena.Examples of atomistic modelling techniques are molecular dynamics (MD), Monte Carlo (MC), and ab initio.Young's modulus for CNTs by utilizing atomistic modelling is reported to be within the range of 0.55 TPa to 1.4 TPa [14-18, 21, 31-35].Atomic modelling is limited to very small scales (length and time) and requires high computing power.Then, continuum modelling methods emerged to overcome these drawbacks.
Continuum modelling assumes that the CNTs have continuous distributions of mass, stiffness, and others.The entire CNT is represented by continuous structures and neglects the lattice structure of the CNT.There are several approaches in continuum modelling of CNTs.A common approach is by utilizing continuum structures in finite element method (FEM) that are coupled with molecular mechanics.Continuum modelling can also be achieved by utilizing analytical methods.Examples of continuum modelling are shell modelling [36][37][38][39][40][41][42], truss modelling [43], spring modelling [19,44], and beam modelling [45].Other approaches include development of 3D continuum model by considering MWCNT as a foliation, dividing of space into continuous stack of leaves [46], tube modelling [47], and hollow cylinder method [48].Young's modulus for CNTs by utilizing continuum modelling is reported to be within the range of 0.4 TPa to 2.52 TPa [19,36,38,39,[41][42][43][44][45][46][47][48].However, care needs to be taken when continuum modelling is utilized, since the lattice structure of the CNT is compromised in this approach.
Nanoscale continuum modelling replaces the carbon to carbon (C-C) bond within the CNT structure with continuum element.Therefore, the lattice structure and the interatomic reactions of CNTs are not neglected in this regard, as opposed to the continuum modelling described above.Properties obtained from atomistic modelling or molecular theories are incorporated into the continuum element which represents the C-C bond.Beam elements are most commonly utilized to replace C-C bond in CNTs [21][22][23][24][25][26][27][28][29][30][49][50][51][52].Spring elements are also utilized by some researchers [20,24,53,54].Other approaches include utilization of continuum shell model coupled with spring constant [55], analytical analysis by utilization of Morse potential, and other nonlinear potentials [25,38,56].Young's modulus for SWCNTs through nanoscale continuum modelling with utilization of beam elements is reported to be within the small range of 0.81 TPa to 1.4 TPa, for wall thickness of 0.34 nm [21][22][23][24][25][26][27][28][29].From the above survey, it is seen that variation of Young's modulus is smallest in nanoscale continuum modelling and this modelling method is more desired than continuum modelling.SWCNTs are formed by rolling graphene sheet into a tube or cylinder.MWCNTs consist of multiple rolled layers (concentric tubes) of graphite sheets.There are three types of configurations for SWCNTs, which are armchair (n, n), zigzag (, 0), and chiral (, ). and  represent the number of unit vectors along two directions in the honeycomb crystal lattice of graphene.These configurations are shown in Figure 1.
Each cell of SWCNTs is formed by C-C bonds, in which each carbon atom is located at prescribed locations, forming hexagon.For 3-dimensional SWCNTs, the thickness of the CNT is considered to be equivalent to the diameter of the carbon atom, which is 0.34 nm.Different thicknesses have been proposed for CNTs as well [30,[49][50][51][52].Since the unit cell of CNT is in the form of hexagon in 2D and polyhedron in 3D, polygonal/polyhedral FEM is found to be ideal in representing the unit cell.Utilization of polygonal FEM in representing CNT can be seen in [57].The author utilized hexagonal elements in 2D and tetrahedral elements in 3D to study heat transfer through CNTs.Work on 3D polyhedral FEM in nanoscale simulation of CNTs is currently not available in the literature.
In this work, nanoscale continuum models of CNTs are developed, by utilizing 3-dimensional virtual node polyhedral elements (VPHEs) [58].The carbon lattice of CNT is represented by 20-node VPHE.Properties of CNT obtained from atomistic simulation are incorporated into the VPHEs.The new CNT models using VPHEs are simulated under tensile load and corresponding Young's modulus is obtained and later compared with those reported in literature.This paper is organized as follows.Structures of the CNTs based on VPHEs are described in Section 2. Simulation setup is shown in Section 3. Results and discussions are provided in Section 4. The paper is finally concluded in Section 5.

New Polyhedral CNTs
Recently, a polyhedral element has been developed, which is termed as VPHE [58].The arbitrary VPHE is partitioned into  number of tetrahedral elements.Each plane of the VPHE is partitioned into triangles which share a common center, .The number of triangles on a face is equal to total number of nodes on the particular face.The fourth node for each triangle (vertex, ) which forms tetrahedral is located at the center of the VPHE.An example of partitioning of an arbitrary VPHE is shown in Figure 2.
One of the advantages of the VPHE is that the element can take arbitrary form with any number of faces, as shown in a patch configuration in Figure 3.The patch consists of six VPHEs in the form of polyhedral elements with 6, 7, and 8 faces.These elements are numbered from 1 to 6 in Figure 3.
A cell of CNT can be represented by the polyhedral element (in nanoscale), by first locating each node of the front plane of VPHE at the center of each carbon within the cell, which forms the hexagonal structure.The overall 3D polyhedral element can then be formed by extruding the front plane by an amount equivalent to the thickness of CNT.Representation of a cell of CNT by the VPHE is shown in Figure 4 (not to scale).20 nodes of the VPHE element are labelled in Figure 4(b).Nodes 1-12 represent the edge nodes and nodes 13-20 represent the common center, with C shared by each triangle on the respective planes/surfaces.
The cell is then arranged in repeating pattern according to specific configuration to form the CNT with certain diameter, length, and thickness.Figure 5 shows an example of CNT with zigzag configuration of (15, 0) formed by the VPHEs.A single VPHE (which is repeated to form the entire CNT structure) is also shown in Figure 5.

Journal of Nanomaterials
A single VPHE is partitioned into 36 tetrahedrons, which are combined together through least square method.CNTs modelled by VPHE elements are governed by following FEM equations (for solid mechanics) [58]: where   V represents stiffness matrix corresponding to tetrahedron  V (3 degrees of freedom per node),  is the Jacobian matrix, , , and  are the natural coordinate system,   V represents partial derivatives of the interpolation functions corresponding to tetrahedron  V ,  represents polyhedral element node,  represents elastic material property matrix,  [V]    represents polyhedral shape function corresponding to tetrahedron  V ,  represents modulus of elasticity, and V represents Poisson's ratio.Once stiffness matrices corresponding to each of the tetrahedrons within a polyhedron (  1 ,   2 , . . .,   V , . . .,    ) are found, the stiffness matrix for a single polyhedron,   , is then calculated by merging the stiffness matrices: Displacements and forces can then be computed for the global system using the following equations: where  represents global force matrix,  represents global stiffness matrix, and  represents global displacement matrix.

Simulation of CNTs
Young's modulus and Poisson's ratio for a single layer graphene are taken as 235.88 N/m and 0.4136, respectively [59].These values are obtained through MD simulations by taking into account the effect of internal lattice relaxation.When an external load is applied, individual atoms of a single cell of CNT undergo displacements which lead to change in the bond lengths and bond angles.These changes lead to unbalanced interatomic forces between the atoms and eventually internal equilibrium is not maintained at each atom.However, the internal equilibrium must be achieved in order to maintain deformed shape of the structure (deformation due to the applied external load) at the macroscopic level.Therefore, individual atoms undergo additional displacements in specific directions internally, which leads to equilibrium without affecting macroscopic deformation.These additional internal displacements of the atoms (or lattices) relax the total potential energy, eventually known as internal lattice relaxation.
Since CNTs are produced by rolling a single graphene sheet, values of Young's modulus and Poisson's ratio can be assigned to the cells of CNT (see Figure 4).For the modelling of CNT through VPHEs, the 2D value of Young's modulus is converted to 3D, by dividing with thickness of the CNT [60], which gives  VPHE = 0.69 TPa.Several FEM simulations are run in the nanoscale by utilizing Mathematica software, for armchair and zigzag configurations.Codes are written to generate coordinates of each node of the CNT cell with armchair and zigzag configurations, with varying length, diameter, and thickness.The length of C-C bond is taken as 0.142 nm.Young's modulus of various CNT models ( CNT ) with different parameters is determined through tensile test and the results are compared with other published data.The tensile test is simulated by fixing one end of the CNT and the other end is subjected to a predetermined displacement of 0.0023 nm (within the elastic region).Net reaction force is calculated by summing all the forces obtained at the fixed end of the CNT (by utilizing (1)).Young's modulus for entire CNT model,  CNT , is then calculated by utilizing following formulas [20,26,28,49,52]: where ∑  represents net reaction force at the fixed end of CNT and  represents prescribed displacement at the free end of CNT. CNT , , , and  represent stiffness, length, diameter, and thickness of the CNT.Effects of CNT diameter , length , and wall thickness  on  CNT are examined, one at a time.Value of the parameter being studied is varied, while the other two parameters are kept constant.

Results and Discussion
Effect of CNT length on  CNT is studied, by varying parameter  for both armchair and zigzag configurations.The length is increased until a steady state solution is obtained.Diameters of the CNT for armchair and zigzag configurations are fixed at 0.8136 nm and 0.46974 nm, respectively, based on  = 6.Thickness of the CNT wall is fixed at 0.34 nm for both configurations.Results of the simulations are shown in Figure 6.It is seen that  CNT increases with length during the initial stage and later converges to a value of 1.23 TPa for armchair and 1.338 TPa for zigzag.The initial variation of  CNT with length (highlighted in Figure 6) is caused by the end conditions.One of the ends of the CNT is fixed and the other end is subjected to load.Therefore, when the length of the CNT is considerably smaller, the boundary condition (at the fixed end) is closer to the applied load and affects the reading.Ratio of length over radius (2/), which is also known as slenderness ratio, should be greater than or equal to 10 to avoid the effect of the end conditions described above [16].Figure 6(a) shows that  CNT starts to  convergetowardsthesteady state reading when  = 10.72 nm, with slenderness ratio of 45.6, while Figure 6(b) shows that the steady state reading is obtained when  = 5.04 nm, with slenderness ratio of 12.4.
Effect of CNT diameter on  CNT is studied by changing the value of coefficient  to 6, 7, and 8, for zigzag and 3, 4, and 5 for armchair configurations, in order to obtain similar diameter between them.Length of CNT is fixed at 6.27 nm for armchair and 10.72 nm for zigzag, with slenderness ratio greater than 10 for all the cases examined.Thickness of CNT is fixed at 0.34 nm.The results are shown in Figure 7. From the results, it is seen that Young's modulus decreases with increasing diameter, for both armchair and zigzag configurations.Armchair configuration yields higher  CNT compared to zigzag configuration. CNT decreases by following the trend of a quadratic curve.This is due to the inversely proportionality relationship of the cross-sectional area of the CNT and  CNT .As observed from (4), larger diameter yields higher cross-sectional area.However, the trend obtained in Figure 7 is in contradiction with some of other employed nanoscale continuum models.Differences in the trends are caused by the type of potential function utilized in the formulations of the interatomic molecular energies.Another factor is the approach taken in modelling and calculation of Young's modulus of the CNTs.These are explained below.
One of the approaches in modelling CNTs in nanoscale is by representing the C-C bond by utilizing continuum element (FEM) such as beam or spring elements.Molecular energies between C-C bonds are represented through interatomic potential functions which are normally utilized in molecular mechanics and molecular dynamics approach.Then, the associated constants are derived by establishing the link between the potential energies and strain energies of the continuum structure [21].There are various types of interatomic potential functions such as Lennard-Jones (LJ), Morse, Brenner, or Tersoff.Potential functions such as Morse and Brenner have been utilized together with spring elements [20,53] and beam elements [28].These potential functions yield increasing value for Young's modulus with increasing CNT diameter [20-23, 26, 28, 49, 53, 54].Closed-loop solution which is developed in [25] by utilizing AMBER force field to determine Young's modulus showed similar increasing trend.However, authors in [61] carried out molecular mechanics simulations by utilizing different potential functions and showed that the trend for Young's modulus versus CNT diameter is influenced by the type of potential function utilized.The authors showed that Young's modulus decreases with increasing CNT diameter when reactive empirical bond order (REBO) potential function is utilized and showed opposite trends for modified Morse potential function and the universal force field (UFF).
Decreasing values of Young's modulus with increasing CNT diameter is obtained when the approach taken in  calculation of Young's modulus of the CNTs is changed as well.For example, in [27], the authors adopted similar approach as in [49] to model SWCNT, but calculation of Young's modulus of the CNTs is modified by taking into account the rigidities in tension and bending.This resulted in decreasing values of Young's modulus with increasing CNT diameter, which is opposite to the results obtained in [49].In [30], the authors calculated Young's modulus of CNTs by comparing the nanoscale continuum model with continuum model (cylindrical solid).The equivalence of energies between these two models was utilized as a framework to develop a new way of calculating Young's modulus.This approach yields decreasing values of Young's modulus with increasing CNT diameter.Similar decreasing trends are also reported in other literatures with different modelling approaches, which takes into account the rigidities in tension or bending of CNTs (such as structural mechanics approach).Such examples can be seen in [32,62] for atomistic modelling, [42,45] for continuum modelling, and [63] for nanoscale continuum modelling.
Effect of CNT wall thickness on  CNT is studied by varying the parameter , for both armchair and zigzag configurations.Length of CNT is fixed at 6.27 nm for armchair and 10.72 nm for zigzag (similar to the previous case).Three sets of CNTs are considered for each configuration, with n = 6, 7, and 8. Results are shown in Figure 8. From the results, it is seen that Young's modulus decreases with increasing wall thickness, for both armchair and zigzag configurations.This is due to the inversely proportional relationship between the cross sectional area of the CNT and  CNT , which is observedfrom (4).Larger thickness yields higher cross-sectional area in (4).Wall thickness of 0.34 nm represents SWCNT, while wall thickness of 0.68 nm represents double walled CNT (DWCNT).These two cases are labelled in Figure 8.It is seen that increasing the number of layers for CNTs do not cause increment in  CNT .
Table 1 summarizes results taken from other references.Results from the current work are in good agreement with other published works.For example, average Young's modulus of 1.2 TPa which is obtained for armchair configuration from this work is close to the results obtained in [3, 14-16, 19, 21].Average value of 1.3 TPa which is obtained for zigzag configuration in current work is similar to the results presented in [16,20].

Conclusions
The VPHEs are shown to be ideal in modelling CNTs with results that match the other published works in literature.This study provides the insight on nanoscale continuum modelling of carbon nanotubes with effect of varying parameters (length, thickness, and diameter) on Young's modulus of the CNTs.Utilization of polyhedral FEM in modelling CNTs with varying parameters is currently not available in literature.Advantage of the method described in this paper is that it is applicable for both 2-dimensional and 3-dimensional analyses.Another advantage is that properties of the VPHEs can be defined accordingly to represent general nanotubes (which contain entities other than carbon).This work is useful in modelling CNTs for composite materials and other future CNT implementations.

Figure 3 :
Figure 3: Patch consisting of six VPHEs in arbitrary form.

Figure 4 :
Figure 4: Representation of a cell of CNT with VPHE: (a) a cell of CNT consisting of 6 carbon atoms, forming a hexagon, and (b) equivalent VPHE element which represents the cell of CNT in 3D.

Figure 7 :
Figure 7: Effect of CNT diameter on Young's modulus of CNT.

Figure 8 :
Figure 8: Effect of CNT wall thickness on  CNT : (a) results for zigzag configuration with different diameters and (b) results for armchair configuration with different diameters.

Table 1 :
Comparative study on Young's modulus.