NONLOCAL CONTINUUM SHELL MODELS FOR TORSION OF SINGLE-WALLED CARBON NANOTUBES by Farzad Khademolhosseini B.Sc., Sharif University of Technology, Iran, 2007 A THESIS SUBMITTED IN PARTIAL FULLFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2009 © Farzad Khademolhosseini, 2009 ii Abstract Carbon nanotubes (CNTs) have attracted much attention from scientists and engineers because of their relevance to a wide range of applications. Various approaches have been used for the characterization of CNT properties, among which continuum modeling has generated much interest due to computational efficiency. However, at the nanoscale the dimensions of a system are comparable to the inter-atomic or inter- molecular spacing of that system, and the material cannot be modeled as a continuum. This is known as the “size-effect”. To overcome the limitations of classical continuum mechanics, modified continuum models have been proposed, among which models based on the concept of nonlocal elasticity have proven effective in quantifying the size- dependent mechanical response of CNTs. This thesis investigates the “small-size” effects in the torsional response of single walled carbon nanotubes (SWCNTs) by developing a modified nonlocal continuum shell model for their torsion. The purpose is to facilitate the design of devices based on CNT torsion by providing a simple, accurate and efficient continuum model that can predict the torsional buckling loads, the frequency of torsional vibrations and the propagation speed of torsional waves. To this end, Eringen’s equations of nonlocal elasticity are incorporated into the classical models for torsion of cylindrical shells given by Timoshenko and Donnell. In contrast to the classical models, the nonlocal model developed here predicts non-dimensional buckling loads that depend on the values of certain geometric parameters of the CNT, allowing for the inclusion of size-effects. In the case of torsional vibrations and propagation of torsional waves, the classical and nonlocal models predict non-dispersive and dispersive behavior, respectively. Molecular dynamics simulations of torsional buckling, axial buckling and torsional vibration of various SWCNTs are also performed, the results of which are compared with the classical and nonlocal models and used to extract consistent values of the nonlocal elasticity constant. Interestingly, the nonlocal elasticity constant depends on the existence of circumferential and/or longitudinal modes in the deformed shape of the CNT. In all loading cases the superiority of the nonlocal model over the classical elasticity model in predicting the size- dependent mechanical response of SWCNTs is established. iii Table of Contents Abstract.............................................................................................................................. ii Table of Contents ............................................................................................................. iii List of Tables ......................................................................................................................v List of Figures....................................................................................................................vi List of Symbols ..................................................................................................................ix Acknowledgements ............................................................................................................x Dedication ..........................................................................................................................xi 1. INTRODUCTION ........................................................................................................1 1.1 Nanotechnology and Carbon Nanotubes .................................................................1 1.2 Carbon Nanotubes as Torsional Elements ...............................................................5 1.3 Methods of Studying Carbon Nanotubes.................................................................7 1.4 Review of Modified Continuum Models for Nanoscale Structures ........................9 1.5 Scope of the Present Work.....................................................................................11 2. NONLOCAL SHELL MODEL FOR BUCKLING OF SWCNTS ........................14 2.1 Review of Eringen’s Nonlocal Elasticity Theory ..................................................14 2.2 Modified Nonlocal Shell Model for Torsional Buckling.......................................16 2.2.1 Modified Timoshenko shell model ..............................................................16 2.2.2 Modified Donnell shell model .....................................................................22 2.3 Modified Nonlocal Shell Model for Axial Buckling .............................................25 2.4 Size-dependency of Nonlocal Models ...................................................................29 3. MOLECULAR DYNAMICS SIMULATION OF CNT BUCKLING...................31 3.1 Review of Molecular Dynamics ............................................................................31 3.1.1 Basic theory .................................................................................................31 3.1.2 Atomic interactions and interatomic potentials ...........................................33 3.1.3 Important factors in molecular dynamics simulations.................................35 3.1.4 Quasi-static molecular dynamic simulations ...............................................37 3.2 Molecular Dynamics Simulation of Torsional Buckling of SWCNTs ..................38 3.2.1 Review of NanoHive ...................................................................................38 3.2.2 Simulation setup ..........................................................................................38 3.2.3 Simulation results ........................................................................................41 iv 3.3 Values of Nonlocal Constant and Shell Thickness................................................46 3.3.1 Review of previous work.............................................................................47 3.3.2 Calculation of e0 and h using optimization..................................................48 3.3.3 Comparison of classical and nonlocal models with MD simulations..........49 3.4 MD simulations of Axial Buckling........................................................................55 3.4.1 Axial buckling simulation setup ..................................................................55 3.4.2 Axial buckling simulation results ................................................................57 3.4.3 Values of nonlocal constant and shell thickness..........................................59 4. NONLOCAL SHELL MODEL FOR TORSIONAL VIBRATION OF CNTS ....61 4.1 Nonlocal Shell Model for CNT Torsional Waves and Vibrations.........................61 4.1.1 Propagating wave solution...........................................................................63 4.1.2 Standing wave / vibration solution ..............................................................65 4.1.3 Comparison of classical and nonlocal models .............................................67 4.2 MD Simulations of Torsional Vibrations ..............................................................71 4.2.1 MD simulation setup....................................................................................71 4.2.2 MD simulation results..................................................................................74 4.3 Comparison of MD Simulation Results with Classical and Nonlocal Models......75 4.3.1 Comparison of radial frequencies ................................................................76 4.3.2 Comparison of group velocities...................................................................78 4.3.3 Comparison of phase velocities ...................................................................82 4.4 Discussion on the Nonlocal Constant Value..........................................................86 5. SUMMARY AND CONCLUSIONS .........................................................................89 5.1 Summary of Present Work and Major Findings ....................................................89 5.2 Suggestions for Future Work.................................................................................91 Bibliography.....................................................................................................................92 Appendix...........................................................................................................................98 v List of Tables Table 3.1 Properties and critical torques and twists (a) zigzag (b) armchair CNTs .... 44 Table 3.2 Values of shell thickness (h) and nonlocal constant (e0) obtained from non- linear least square fitting of MD simulation results with classical shell and nonlocal shell models................................................................................... 49 Table 3.3 Properties and critical forces and strains of several armchair CNTs ........... 58 Table 3.4 Values of shell thickness (h) and nonlocal constant (e0) obtained from non- linear least square fitting of MD simulation results for axial buckling ....... 59 Table 4.1 Vibration periods and natural frequencies of several different mode-shapes of (6,6) armchair CNT ................................................................................. 74 Table 4.2 Vibration periods and natural frequencies of several different mode-shapes of (10,10) armchair CNT ............................................................................. 75 Table 4.3 Surface shear modulus, surface density and shear wave speed c for (6,6) and (10,10) armchair CNTs ........................................................................ 76 Table 4.4 Comparison of group velocities of classical and nonlocal models with MD simulation results for (a) (6,6) and (b) (10,10) armchair CNTs................... 79 Table 4.5 (a) Comparison of phase velocities of classical and nonlocal models with MD simulation results for (6,6) armchair CNT ........................................... 82 (b) Comparison of phase velocities of classical and nonlocal models with MD simulation results for (10,10) armchair CNT ....................................... 83 Table 4.6 Obtained values of nonlocal constant .......................................................... 88 vi List of Figures Figure 1.1 Definition of chiral indices (m,n), unit vectors a1 and a2, rollup vector Ch and chiral angle α........................................................................................... 2 Figure 1.2 Structural configurations of zigzag, armchair and chiral nanotubes.............. 3 Figure 1.3 Cross section of single, double and triple-walled carbon nanotubes ............. 3 Figure 1.4 Classification of different carbon nanotube structures .................................. 4 Figure 1.5 Schematic of an electromechanical paddle oscillator with SWCNT spring .. 6 Figure 1.6 Various methods used to study the mechanical properties of carbon nanotubes ....................................................................................................... 8 Figure 2.1 Cylindrical shell representation of a SWCNT with the coordinate system used .............................................................................................................. 16 Figure 2.2 Definition of forces and moments per unit length acting on a shell ............ 18 Figure 2.3 Definition of classical and nonlocal buckling loads for different values of the nonlocal constant e0 (λ=0.5 and n=2 have been used) ........................... 30 Figure 3.1 Summary of steps in the molecular dynamics simulation technique ........... 32 Figure 3.2 Definition of bond length and bond angles used in molecular mechanics... 34 Figure 3.3 A uniform counterclockwise twist of 1 radian induced in a (5,5) armchair CNT by mathematical rotation of atomic coordinates around the central axis (Atoms marked with black circles are fixed in space) ................................. 39 Figure 3.4 Relaxation of a (5,5) armchair nanotube with end twist of 3 radians corresponding to a uniform shear strain of 0.1 (a) Start of simulation (b) After relaxation for 6000 time-steps (c) After relaxation for 6500 time-steps (d) After relaxation for 40000 time-steps (e) Progression of potential energy of the CNT during the 40000 timesteps of relaxation.................................. 40 Figure 3.5 (a) Representation of predicted buckling mode-shape using MATLAB (b) vii Actual buckling mode-shape predicted using MD simulations ................... 42 Figure 3.6 MD simulation results (a) Variation of potential energy of (8,0) nanotube with applied shear strain (twist) (b) Stress – strain curve calculated using slope of potential energy curve .................................................................... 43 Figure 3.7 Comparison of buckling torques from classical shell models with MD simulation results for (a) zigzag and (b) armchair CNTs based on properties given in Table 3.2......................................................................................... 50 Figure 3.8 Comparison of buckling torques from classical and nonlocal Timoshenko shell models with MD results for (a) zigzag CNTs and (b) armchair CNTs based on properties given in Table 3.2 ........................................................ 52 Figure 3.9 Comparison of buckling torques from classical and nonlocal Donnell shell model with MD results for (a) zigzag CNTs and (b) armchair CNTs based on properties given in Table 3.2................................................................... 53 Figure 3.10 A compressive axial strain of 0.04 is induced in a (10,10) armchair CNT by manually changing the coordinates of carbon atoms (atoms at both ends of CNT shown in dark color are fixed) ............................................................ 56 Figure 3.11 Different instances of relaxation of a (10,10) armchair nanotube under an axial compressive strain of 0.04. (a) simulation start (b) after relaxation for 19250 time-steps (c) after relaxation for 20250 time-steps (b) after relaxation for 30000 time-steps.................................................................... 56 Figure 3.12 Axial buckling mode-shape of (10,10) armchair nanotube at critical strain 57 Figure 3.13 Progression of potential energy with axial buckling strain of a (10,10) armchair nanotube........................................................................................ 58 Figure 3.14 Comparison of axial buckling loads from (a) classical and (b) nonlocal shell model with MD results for several armchair CNTs based on properties given in Table 3.3 .................................................................................................. 60 Figure 4.1 Comparison of classical and nonlocal models at different wavenumbers and viii for different values of the nonlocal constant e0 (a) phase velocities (b) group velocities ...................................................................................................... 70 Figure 4.2 Free torsional vibrations of a (6,6) armchair CNT at the eighth natural modes-shape corresponding to k′=8, where the mode-shape is induced by mathematically changing the atomic coordinates ........................................ 72 Figure 4.3 Free torsional vibrations of a (6,6) armchair CNT at the 8th natural modes- shape corresponding to k′=8. . Different instances of half a vibration period (τ /2=160 fs) are shown................................................................................ 73 Figure 4.4 Comparison of local and nonlocal models for the normalized radial frequency versus mode-shape number with MD simulation results (a) (6,6) CNT (b) (10,10) CNT .................................................................................. 77 Figure 4.5 Comparison of local and nonlocal models for the normalized group velocity versus wavenumber with MD simulation results (a) (6,6) CNT (b) (10,10) CNT.............................................................................................................. 81 Figure 4.6 Comparison of local and nonlocal models for the normalized phase velocity versus wavenumber with MD simulation results (a) (6,6) CNT (b) (10,10) CNT.............................................................................................................. 85 ix List of Symbols a Mean radius of cylindrical shell used to model a CNT d Internal characteristic length (C-C bond length) e0 Nonlocal constant NT crF Nonlocal Timoshenko axial buckling load of cylindrical shell T crF Classical Timoshenko axial buckling load of cylindrical shell h Thickness of cylindrical shell used to model a CNT l Length of cylindrical shell used to model a CNT NT crM Nonlocal Timoshenko buckling torque of cylindrical shell T crM Classical Timoshenko buckling torque of cylindrical shell ND crM Nonlocal Donnell buckling torque of cylindrical shell D crM Classical Donnell buckling torque of cylindrical shell m Longitudinal half-wave number n Circumferential wave number u Displacements in the axial / longitudinal direction of cylindrical shell v Displacements in the circumferential direction of cylindrical shell w Displacements in the radial direction of cylindrical shell ξ Nonlocal parameter ξ=(e0.d)2 λ Longitudinal wave length λ=mπa/l υ Poisson’s ratio ( 0.2υ  for CNTs) x Acknowledgements I am heartily thankful to my supervisor, Professor Nimal Rajapakse, and my co- supervisor, Dr. Alireza Nojeh, whose encouragement, guidance and support from the initial to the final level enabled me to develop a deep understanding of the subject. I would like to show my gratitude to Dr. Srikanth Phani, for fruitful discussions and for his valuable inputs to this thesis. I would also like to express my gratitude to Professor Farrokh Sassani for being a member of the defense committee and for his useful comments on the thesis. I gratefully acknowledge the support from Professor Nimal Rajapakse’s grant from the Natural Sciences and Engineering Research Council of Canada (NSERC). Last but no least, I am indebted to my family, specially my parents, for their undying support and constant encouragement. xi gÉ Åç ÑtÜxÇàá 1 Chapter 1 INTRODUCTION 1.1 Nanotechnology and Carbon Nanotubes Nanotechnology [1] is a fast growing field of science and technology which deals with the development of materials, structures and devices at the nanoscale. A nanometer is one billionth of a meter. Nanoscale materials are those with at least one dimension in the nanometer range (between 1 to 100 nanometers). At these extremely small dimensions, quantum effects are significant and the discrete nature of matter gives rise to interesting phenomena which are rarely observed at the Macro-scale. These interesting properties can be controlled, manipulated and used in the development of novel devices [2, 3]. One interesting nanostructure that has attracted a lot of attention is the carbon nanotube (CNT). The discovery of CNTs in the 1990s [4] is regarded as a revolutionary step in advancement of nanotechnology. Based on their unique structure and large aspect ratios, different properties were hypothesized for CNTs, some of which have been proven and are currently utilized in high-tech applications. A fundamental understanding of CNTs, their unique structure and exceptional characteristics is beneficial to advancement and innovation in the frontiers of nanoscience and technology. So what exactly are carbon nanotubes? Just like diamond and graphite, a carbon nanotube is an allotrope of carbon. It is a single molecule of carbon atoms which can be thought of as a graphene sheet which has been rolled up into a cylindrical shape. Note that this is not the way nanotubes are made in practice; however, the above description helps illustrate the nature of the CNT structure. Carbon nanotubes can have length to diameter ratios of up to 28,000,000:1 [5], which is significantly higher than any other material known. Such high aspect ratios at nanoscale dimensions essentially cause carbon nanotubes to behave as one-dimensional structures and this contributes to their outstanding properties. Based on the direction in which a graphene sheet is rolled up to form a carbon nanotube, different CNT structures are possible and so there is a need to differentiate 2 between them. This is achieved by using Chiral indices (m,n) to name different CNT structures [6]. Defining an orthogonal coordinate system with unit vectors a1 and a2 in the plane of a graphene sheet, the rollup vector Ch equals to ma1+na2 (Figure 1.1). Once the graphene sheet has been rolled up into the cylindrical form, the roll-up vector Ch constitutes the perimeter of the carbon nanotube. Furthermore, the CNT axis is perpendicular to the direction of the roll-up vector. The angle α between the roll-up vector Ch and unit vector a1 is called the “Chiral angle”. CNTs with a chiral angle equal to zero (or n=0) are called ‘Zigzag’, those with a chiral angle of 30º (or n=m) are called ‘Armchair’, and those with chiral angles between 0º and 30º (or 01 are considered. Substitution of equations (2.37) into the modified governing differential equations (2.36) and solution of the resulting eigenvalue problem yields the following relation for the modified critical axial buckling load of a cylindrical shell: 2 2 22 2 2 2 2 0 2 2 1. 13(1 ) .2 1 ( ) NT cr cr Eh n n F ah de n a π υσ π λ − +−= = + + , (2.38) where λ=mπa/l. The non-dimensional form of the critical buckling load is: 2 22 2 2 2 2 2 0 2 2 1. 13(1 ) 1 ( ) NT NT cr cr n nFF Eh de n a π υ λ − +−= = + + . (2.39) For comparison the classical Timoshenko relations for the critical axial buckling load and its non-dimensional form are: 28 2 2 22 2 1. 13(1 ) T T cr cr F nF Eh n π υ −= = +− . (2.40) Therefore, 2 2 2 2 0 21 ( ) T T cr cr NT NT cr cr F F de n F F a λ= = + + . (2.41) Comparison of equations (2.39) and (2.40) clearly shows that the former is size- dependent. This size-dependency is further explained in the next section. Equations (2.38), (2.39) and (2.40) are based on the critical axial buckling wavelength calculated as [43]: 2 2 2 2 ( 1) 2 3(1 ) hn n a λ υ −= − (2.42) 29 2.4 Size-dependency of Nonlocal Models Looking back at the classical relations for buckling loads derived in equations (2.18), (2.31) and (2.40) it is easily seen that the non-dimensional form of these buckling loads are all independent of size. In other words, regardless of the geometric dimensions and aspect ratios of a cylindrical shell, classical models always predict the same value for the non-dimensional buckling load. In contrast, looking at the nonlocal forms of buckling loads derived in equations (2.17), (2.30) and (2.39) it is easily seen that the non- dimensional form of these buckling loads are all size-dependent. This is due to the presence of geometric ratios in the denominators of the aforementioned equations. As an example the denominator of equation (2.39) is analyzed here: 2 22 2 2 2 2 2 2 2 0 021 ( ) 1 d d de n e n m a la λ π⎛ ⎞⎛ ⎞ ⎛ ⎞+ + = + +⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎝ ⎠ ⎝ ⎠⎝ ⎠ . (2.43) The magnitude of the expression of equation (2.43) depends on the value of two geometric ratios, d/a and d/l, where d is the inter-atomic distance or the C-C bond length and a and l are the radius and length of the CNT, respectively. Since d is constant or almost constant for all CNTs (d ≈ 1.41 angstroms), based on their respective length and diameter size, different CNTs will have different d/a and d/l ratios and thus different non- dimensional buckling loads. As the lengths and diameters of CNTs become smaller and comparable to the inter-atomic distance d, the d/a and d/l ratios become larger and have a greater impact on the magnitude of equation (2.43). This is how nonlocal models are able to account for “small” size effects. The magnitude of the expression of equation (2.43) also depends on: i) The buckling mode-shape through the values of m (longitudinal half-wave number) and n (circumferential wave-number). As wave-numbers increase, the corresponding wave-lengths decrease and become comparable to the inter- atomic distance d. ii) The value of the nonlocal constant e0, which has to be determined for each material independently. 30 Figure 2.3 which is based on equations (2.32) and (2.41) gives a comparison of classical and nonlocal buckling loads for different values of the nonlocal constant e0. 2 4 6 8 10 12 14 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 CNT Radius (Angstrom) C la ss ic al B uc kl in g Lo ad / N on lo ca l B uc kl in g Lo ad e 0 =0.8 e 0 =0.6 e 0 =0.4 e 0 =0 Figure 2.3- Comparison of classical and nonlocal buckling loads for different values of the nonlocal constant e0 (λ=0.5 and n=2 have been used). It is seen that the value of the nonlocal constant can significantly affect the buckling loads predicted by the nonlocal models. The difference between classical and nonlocal models is significant for CNTs with small diameters, but this difference becomes negligible (less than 5%) at larger diameters (above 3 angstroms) where the nonlocal and classical models converge. Based on the above observations it is important to know the buckling mode-shape (values of m and n) and the correct value of the nonlocal constant e0 in order to calculate the magnitude of small-size (nonlocal) effects on buckling loads. This task is undertaken in the next chapter where molecular dynamics (MD) simulations will be used in conjunction with classical and nonlocal models to calculate the values of the aforementioned parameters. 31 Chapter 3 MOLECULAR DYNAMICS SIMULATION OF CNT BUCKLING 3.1 Review of Molecular Dynamics Molecular dynamics is an atomistic/molecular simulation technique which was developed in the late 50s [45]. In recent years with the fast growth of emerging fields such as nanotechnology, molecular dynamics has provided researchers with a simple yet accurate method for the study of atomic and molecular systems. Molecular dynamics calculates the motions of atoms or molecules in a system. The time evolution of interacting atoms is followed by integrating their equations of motion which follow the laws of classical mechanics, most notably Newton’s second law of motion, F=ma [46]. Due to the important role of MD simulations in this work and its importance to nanoscale research in general, the basic concepts behind the molecular dynamics simulation method are presented here. 3.1.1 Basic theory The classical molecular dynamics simulation technique follows a simple and systematic approach: i. For a system or set of atoms and molecules, an original spatial distribution is assumed. ii. Based on this spatial distribution a potential field is defined for this set of atoms using empirical, semi-empirical or first-principles approaches. iii. Knowing the spatial coordinate of each particle in the system the spatial gradient of the potential field is calculated and utilized in finding the force acting on every particle in the system. iv. Knowing the mass of each particle in the system and the force exerted on it, Newton’s second law is used to calculate the acceleration of each and every particle. 32 v. All particles are then be allowed to move with their respective accelerations for a specific period of time (timestep) and their respective displacements during that timestep are calculated. vi. These displacements are used along with the previous spatial distribution to find the new spatial distribution of the atoms in the system. vii. Steps (ii) to (vi) are repeated to simulate the evolution of the system for a specific amount of time or until a desired state (such as the equilibrium configuration of the system) is reached. Figure 3.1 shows a block diagram with a summary of steps in a MD simulation as described above. Figure 3.1- Summary of steps in the molecular dynamics simulation technique. Assume an initial spatial distribution for the particles in the system Define the potential field created by the particles Calculate forces acting on particles Calculate the acceleration of each particle Calculate displacements of particles for one timestep due to their respective accelerations Calculate new spatial distribution of particles in the system 33 3.1.2 Atomic interactions and inter-atomic potentials As previously stated in a molecular dynamics simulation the forces exerted on the particles are derived from a potential function which depends on the particle coordinates. For a system of n particles the force exerted on the ith particle is [46]: 1 2( , ,..., )i i nF V r r r= −∇ , (3.1) where V is the potential function; a description of the terms by which the particles in the simulation interact. The most common classical potentials are based on molecular mechanics [47], which model the particle-particle interactions due to structural and conformational changes but are unable to reproduce chemical reactions. Classical potentials are based on the following two simplifying assumptions: - Due to their comparatively small mass and fast dynamics electrons respond instantaneously to the motion of their nuclei. - Nuclei can be considered as point particles (point masses) which follow the laws of classical Newtonian dynamics. In molecular mechanics the potential energy of a system of atoms is the sum of covalent and non-covalent bonding energies [47]. co alent non co alentE E Eν ν−= + (3.2) The covalent bonding energies are a function of the bond length (pair-interaction) and different bond angles (many-body interactions), specifically the in-plane, out of plane and dihedral angles. The non-covalent bonding energies are the sum of potential energies due to electrostatic and van der waals interactions; co alent bond in plane angle out of plane angle dihedral angleE E E E Eν − − −= + + + (3.3) non co alent electrostatic vanderwaalsE E Eν− = + (3.4) 34 Figure 3.2 summarizes the different distances and angles which affect the potential energy through the interactions of two or more atoms. r α the bond length r in-plane angle α φ θ Out of plane angle φ Dihedral angle θ Figure 3.2- Definition of bond length and bond angles used in molecular mechanics If a finer level of detail is required for the MD simulation, potentials based on quantum mechanics can be used. Although quantum mechanical potentials produce high levels of accuracy they are often complicated and computationally expensive. Currently there are many different semi-empirical, empirical and quantum mechanical potentials available for use in MD simulations. Some are pair potentials like the non-bonded Lennard-Jones potential [48] used for calculating the van der waals forces. Others are many body potentials like the Tersoff [49] potential originally used in simulating systems consisting of carbon and semi-conductor atoms. Also, potential functions have been developed which allow the creation and termination of atomic bonds by incorporating effects of bond-order [50]. The correct choice of what potential to use in a MD simulation depends on the nature of the simulation, the material or materials being simulated and the trade-off between accuracy and computational efficiency. 35 3.1.3 Important factors in molecular dynamics simulations Just like any other experimental or numerical modeling technique, the accuracy and efficiency of molecular dynamics simulations depend on several factors, which have to be addressed carefully. In MD simulations used to characterize the mechanical properties of systems, factors such as temperature control, potential function, timestep and strain rate are very important and should be addressed with great care. Temperature control Most molecular dynamics simulations are used to visualize the behavior of atomistic systems under specific temperature constraints and so there is a need to control the temperature of the system being simulated. It has been shown that variations in temperature affect the mechanical properties of nanostructures including carbon nanotubes [51]. Temperature control in MD simulations is done through thermostats. Temperature is a measure of the average kinetic energy of the system. A thermostat controls the temperature of a system by controlling its average kinetic energy. For example, in a constant temperature simulation, after each step (or a specific number of steps) of a MD simulation, the velocities of a group of atoms or all the atoms in the system are scaled in such a way that the average kinetic energy remains constant. Different thermostatic functions can be used in MD simulations based on the accuracy and efficiency that is needed. One famous thermostat used in MD is the Brendson thermostat [52], which is explained in detail in section 3.2.1. Timestep The timestep is the length of time between two consecutive steps in a MD simulation. The value chosen for a timestep is of significant importance to the results and functionality of a MD simulation. The timestep controls the trade-off between accuracy and computational efficiency in a MD simulation. Larger timesteps increase the computational efficiency and smaller timesteps increase the accuracy of the simulation. If a timestep is chosen too large, the system might become unstable. As a general rule of thumb, the timestep chosen for a MD simulation should be smaller than one-tenth of the vibration period of an atom or particle in that system [53]. 36 Potential function and potential field One of the most important factors affecting the accuracy of a molecular dynamics simulation is the accuracy and completeness of the potential function/field being used to describe the interaction of particles in the system [53]. The accuracy of a potential function can be increased in many ways, for example by including higher order terms which account for many-body interactions, or by accounting for long-range forces and non-bonded interactions between particles. However, as potential functions become more complicated, the computational time of a simulation also increases and thus there is a trade off between the accuracy and computational efficiency of a simulation. It is, therefore, important to analyze a system, see which inter-particle interactions could have a significant effect on its behavior and choose a potential function which addresses those interactions. Strain rate The time rate of displacement or the strain rate of an imposed deformation can affect the mechanical properties obtained from molecular dynamic simulations of nanostructures. Properties such as the yield or failure point and yield or failure mechanism are considerably dependent on the strain rate of a molecular dynamics simulation [51, 54]. At lower strain rates the system has more time to relax to an equilibrium state and it is more probable that the system deforms along the minimum energy path, thus producing more accurate results. It has been shown that at lower strain rates carbon nanotubes yield at lower values of strain [51]. The failure mechanism of CNTs is also affected by the strain rate [51]. Generally simulations performed at lower strains give more accurate results by allowing more time for the interaction of the particles in the system. Practical strain rates are impossible to simulate using conventional MD. As an example, let us look at the MD simulation of a constant strain tension test on a CNT. To induce a 10% strain in the undeformed CNT, a typical MD simulation with 100000 timesteps of 1 fs each, would simulate an extremely high strain rate of 109 (1/s). 37 3.1.4 Quasi-static molecular dynamic simulations One way of addressing the strain rate effect on the mechanical properties of nanostructures calculated using molecular dynamics is the use of quasi-static simulations. In contrast to conventional molecular dynamics simulations which are strain rate dependent, quasi-static simulations are independent of strain rate. They assume that mechanical deformations take place simultaneously throughout the nanostructure, simulating conditions of infinitely small strain rates. This is better understood in the following example. In a conventional molecular dynamics simulation of CNTs undergoing axial tension, the strain is applied by moving the carbon atoms at one end of the CNT with constant speed and allowing the rest of the atoms in the system to follow due to inter- atomic interactions. This results in a constant strain rate simulation; however, due to the very small timesteps required for a MD simulation these strain rates are very large. As discussed in section 3.1.3, simulations performed at lower strain rates are more accurate and an ideal simulation is one that can simulate very small strain rates. The strain rates in a conventional molecular dynamics simulation can be reduced by reducing the speed of the imposed displacement but this is limited by the time required to simulate a certain amount of displacement. That is, reducing the strain rate increases the time and the number of time steps required to achieve a certain strain and thus decreases the computational efficiency of a conventional molecular dynamics simulation. In quasi- static molecular dynamics simulation of a CNT under axial tension, the axial deformation is induced on the system by mathematically changing the spatial position of carbon atoms to an assumed and pre-defined axial deformation and then allowing the system of atoms to relax to their equilibrium position under the imposed axial strain. In this way, the quasi-static molecular dynamics may overcome the strain rate effect and helps simulate deformations taking place at infinitely small strain rates. However, quasi-static simulations need an assumption on the general deformed shape of the system and thus should only be used for simple simulations where the general deformation can be predicted based on the loading conditions. 38 3.2 Molecular Dynamics Simulation of Torsional Buckling of SWCNTs In this section the simulations performed on the torsional buckling of different types of single-walled carbon nanotubes are presented. First, a brief review of the MD simulator and the simulation setup is given and then the simulation results are presented and compared to previous numerical and experimental results where available. 3.2.1 Review of NanoHive The simulations presented in this work were all performed using the molecular dynamics simulator “NanoHive”. NanoHive is a free open source MD simulator which has certain features that can be used to model different loading conditions of nanostructures including CNTs; however, it has some limitations which will be discussed later in the related sections. NanoHive requires three input files for each simulation. A “nml” file which contains the type and spatial coordinates of all atoms in the system, and where different atoms can be defined as part of a group using the “atomset” definition. A “xml” file which contains information such as the timestep, temperature and duration of the simulation, other information such as the thermostat and potential functions used, as well as the loads or velocities which should be applied to different atomsets. Finally, a “tcl” file showing the simulation flow and defining what outputs are required. NanoHive produces different data as output to a simulation. A text file is produced which contains the potential, kinetic and total energies of the system as well as the temperature at selected timesteps during the simulation. It is also possible to output text files containing the spatial coordinates of the atoms in the system at any desired timestep. Another output of the NanoHive simulator is a “nc” file which contains time histories of the displacement of the atoms in the system. This file can be used as the input to the NanoHive post-processor called “Hivekeeper”, which provides a visual display of the evolution of the system with time. 3.2.2 Simulation setup Quasi-static molecular dynamics simulations are performed on armchair and zigzag nanotubes of different diameters. The process used to model the torsional buckling 39 and find the buckling load is presented here. To determine the buckling state of a certain nanotube under torsion, the modeling starts from an initial non-twisted configuration followed by uniform twists (φ) applied incrementally along the nanotube (Figure 3.3). The new atomic coordinates are subsequently used as the input to the molecular dynamics simulator (NanoHive). The atoms at both ends of the twisted CNT are then fixed in space (Figure 3.3) to simulate the pined-pined boundary condition and MD is used to perform relaxation on the twisted CNT until equilibrium is reached and the potential energy of the system converges to a minimum value. Figure 3.3- A uniform counterclockwise twist of 1 radian induced in a (5,5) armchair CNT by mathematical rotation of atomic coordinates around the central axis (Atoms marked with black circles are fixed in space). Based on the discussion presented in section 3.1.3 [53] a timestep of 0.5 fs is chosen for the molecular dynamics simulations presented here. By looking at the evolution of the system energy with respect to time, it is found that for the nanotubes simulated here, 40000 timesteps are enough to reach equilibrium. For each nanotube, the above simulation is conducted at different values of twist / shear strain and it is seen that above a certain value of twist, which is identified as the critical twist for buckling (φcr), the nanotube collapses into a buckled shape when allowed to relax for a sufficient amount of time (Figure 3.4). Mathematical rotation of atomic coordinates 40 Figure 3.4- Relaxation of a (5,5) armchair nanotube with end twist of 3 radians corresponding to a uniform shear strain of 0.1 for 40000 timesteps (a) Start of simulation (b) After relaxation for 6000 timesteps (c) After relaxation for 6500 timesteps (d) After relaxation for 40000 timesteps (e) Progression of potential energy of the CNT during the 40000 timesteps of relaxation All simulations are performed using the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential [55]. The AIREBO potential is an extension of the commonly used REBO potential developed for solid carbon and hydrocarbon molecules [50]. The extensions that AIREBO potential provides are: a) Non-bonded, intermolecular interactions, modeled with a Lennard-Jones (LJ) 12-6 potential 41 b) Four-body torsional interaction, modeled with a torsional potential with a single minimum A detailed description of the AIREBO potential is presented in appendix A. In the present simulations, once the buckling has occurred, the four-body torsional potential (which includes the dihedral angle change due to bond torsion) becomes important and can affect the equilibrium deformation of the system. The torsional buckling simulations were all performed at room temperature, i.e., 300 Kelvin. NanoHive uses the van Gunstern-Berendsen thermostat [52] which has the following form: 1 2 01 11( ) 2 T Tt T t t λ τ ⎡ ⎤⎛ ⎞⎢ ⎥⎜ ⎟Δ= + −⎢ ⎥⎜ ⎟⎢ ⎥⎜ ⎟− Δ⎢ ⎥⎝ ⎠⎣ ⎦ (3.5) where λ is the temperature scaling factor, Δt is the timestep, τT is a time constant and T0 is the ambient temperature. The scaling factor is used after each iteration to scale the velocities and maintain a near constant temperature around the target temperature T0. 3.2.3 Simulation results Before comparing the results for the critical torques from MD simulations with the nonlocal elasticity shell models, the buckling mode-shapes obtained from the MD simulation are compared with the assumed mode-shapes of the continuum models (e.g. equation (2.16)). For a (8,0) zigzag nanotube with a length of 10 nm, a radius of 0.31 nm, a wall thickness of 0.08 nm and a Poisson’s ratio of 0.2 [56], the buckling wavelength calculated from equation (2.20) corresponds to λ=0.75→ m=7.4 & m/2=3.7. This means that there are 3.7 wavelengths in the longitudinal direction of the buckled nanotube. Figure 3.5 shows a MATLAB representation of the assumed buckling mode-shape of the above nanotube based on equations (2.16) and (2.20) and the actual mode-shape obtained from the MD simulations. 42 Figure 3.5- (a) Representation of predicted buckling mode-shape using MATLAB; (b) Actual buckling mode-shape predicted using MD simulations. Taking into consideration the effects of the fixed boundary conditions at the nanotube ends in the MD model, it is clear from Figure 3.5 that buckling mode-shapes of the nonlocal continuum and MD models agree closely. The MD simulator allows the calculation of the potential, kinetic and total energies of the nanotube, its temperature and coordinates of the atoms at each time step. Figure 3.6 (a) shows the potential energy (U) progression of a (10,0) nanotube based on 18 simulations corresponding to different twist angles / shear strains. The slope of the potential energy – twist curve can be related to the torque to obtain a torque-twist relationship. The following basic relationships exist between U, φ, M and torsional stiffness (K): 2 2 ( )( ) ( )( ) dUM d d UK d ϕϕ ϕ ϕϕ ϕ = = (3.6) For a thin cylindrical shell of length L, wall thickness h and shear modulus G the torsional stiffness is defined by 32K a Gh Lπ= . Therefore: 2 3 2 1 ( ). 2 d UG h a d ϕ π ϕ= (3.7) 43 0 2 4 6 8 -1.35 -1.345 -1.34 -1.335 x 10-15 Shear Strain (%) To ta l P ot en tia l E ne rg y (J ) Figure 3.6- MD simulation results a) Variation of potential energy of (8,0) nanotube with applied shear strain (twist), b) Torque vs strain curve calculated using slope of potential energy curve. Equation (3.7) can be used to calculate the surface shear modulus (G.h) directly from the results of MD simulations and without the need for any assumptions on the value of wall thickness h. (a) (b) 44 It is important to note that based on the length to diameter aspect ratio (L/d), a CNT may buckle in the shaft-shape or shell-shape pattern [56]. Previously, the classical formula of equation (2.18) has mostly been used to study the shell-shape buckling of CNTs, which occurs when [56], 0.5 0.2532 2 3.9 (1 ) 1 4 3 a L D h d Eha υ π −⎛ ⎞ ⎛ ⎞− ⎛ ⎞⎜ ⎟ ≤ ≤ ⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠⎝ ⎠ (3.8) A series of torsional buckling simulations are performed on zigzag and armchair nanotubes with different diameters using the MD software. The lengths of the simulated nanotubes are chosen to keep the aspect ratio approximately constant and within the shell-shape buckling range of equation (3.8). Table 3.1 shows the critical buckling torque (Mcr) and critical buckling twist (φcr) of several armchair and zigzag CNTs. Table 3.1- Properties and critical torques and twists of: (a) zigzag (b) armchair CNTs. Chiral Indices Diam.(Ǻ) Length (Ǻ) L / Diam G.h(GPa nm) Mcr (N.m) φ cr (radians) (10,0) 7.75 122.0 15.7 136.15 6.53E-18 1.65 (12,0) 9.30 140.1 15.1 131.19 7.22E-18 1.22 (14,0) 10.8 169.8 15.2 127.28 7.67E-18 1.05 (16,0) 12.4 184.7 14.9 123.34 8.20E-18 0.85 (20,0) 15.5 245.8 15.9 116.00 9.00E-18 0.71 (a) Chiral Indices Diam.(Ǻ) Length (Ǻ) L / Diam G.h (GPa nm) Mcr (N.m) φ cr(radians) (5,5) 6.70 98.91 14.8 159.11 5.70E-18 1.50 (6,6) 8.05 118.4 14.8 153.17 6.45E-18 1.24 (7,7) 9.40 139.9 14.9 150.70 7.38E-18 1.05 (8,8) 10.7 159.9 14.9 148.80 8.01E-18 0.90 (8,8) 10.7 243.0 22.7 148.50 8.04e-18 1.38 (8,8) 10.7 319.3 29.9 148.60 8.02e-18 1.81 (10,10) 13.4 199.0 14.8 145.32 9.66E-18 0.70 (12,12) 16.1 238.1 14.8 144.12 1.09E-17 0.55 (b) 45 It is evident from Table 3.1 that armchair nanotubes generally have a higher torsional stiffness compared to zigzag nanotubes of the same dimensions (this trend has been reported previously [57]). However, armchair nanotubes buckle at a lower strain and their critical buckling torque is only slightly higher than zigzag nanotubes of similar diameters. The shear modulus values obtained from the present MD simulations are within the range of previously reported experimental and numerical results [6, 28]. Using MD simulations, Yakobson et al. [58] reported a value of 363 GPa.nm for the in-plane stiffness (E.h) and a Poisson’s ratio of 0.19 for a (7,7) armchair CNT. Using these two values the surface shear modulus can be calculated as G.h=E.h/2(1+ν)=152.5 GPa.nm. The surface shear modulus obtained in the present work for a (7,7) armchair nanotube (Table 3.1(b)), is 150.7 GPa.nm. This is almost identical to the value reported by Yakobson et al. (less than 1.5% difference). The buckling strains derived here are lower than those reported by Yakobson et al. [58]. This could be due to the following three reasons; - The MD simulations presented here use the AIREBO potential instead of the Tersoff-Brenner potential used by Yakobson. The AIREBO potential is an improvement to the Tersoff-Brenner potential. - It is well known that MD simulation results are strain-rate dependent. As the strain rate decreases and the system has more time to relax to an equilibrium state, it is seen that properties such as the yield strain decrease [51, 54]. In this work quasi-static simulations are performed which consider deformations taking place at infinitely small strain rates. Quasi-static simulations make sure that the deformation of the tube happens simultaneously at all points and there are no localization effects. This contributes to the fact that our buckling strains are smaller than those reported by Yakobson. - The simulation timestep has a significant influence on the buckling strains calculated using a MD simulation. The timestep chosen has to obey the 10% rule discussed in section 3.1.3, which means for a CNT the timestep has to be smaller than 0.8 fs. This fact is verified through the MD simulations carried out in this work, as it was seen that increasing the simulation timestep from 0.5 fs to 1 fs 46 increased the obtained buckling strain of a (5,5) armchair nanotube from 0.051 to 0.067. That is an increase of more than 30 percent which could lead to a large amount of error. In comparison, reduction of the timestep from 0.5 fs to 0.25fs did not have any considerable effects on the simulation results, as the buckling strain calculated using the 0.25 fs timestep was 0.051, the same as what is seen for the 0.5fs timestep. The higher values of buckling strain derived by Yakobson et al. could be due to using larger values for the timestep. Hall et al. have performed experimental measurements on the torsional properties of single-walled CNTs [13] . Assuming a wall thickness of 3.4 Angstroms, they calculated an average shear modulus of 0.455 TPa for single-walled carbon nanotubes with an average diameter of 1 nm. This corresponds to a surface shear modulus (G.h) of 154.7 GPa.nm. Results presented in Table 3.1(b) demonstrate calculated values of about 150 GPa.nm for the surface shear moduli of armchair CNTs with diameters around 1 nm. This suggests that the CNTs used by Hall et al. are predominantly comprised of armchair SWCNTs and once again, the accuracy of the simulations presented in this work is verified. 3.3 Values of the Nonlocal Constant and Shell Thickness An important issue in the application of the nonlocal elasticity models to CNTs is the determination of the values of the nonlocal elasticity constant e0, equivalent shear modulus and shell thickness. Use of the modified nonlocal models derived in equations (2.17) and (2.30) is only possible if the correct values of parameters such as the surface shear and surface Young’s modulus (G.h and E.h), buckling wave-numbers m and n, Poisson’s ratio ν, nonlocal constant e0 and shell thickness h are known. As seen in section 3.2.3, the correct values for G.h, E.h, m and n can be defined in terms of measurable characteristics of a CNT without any assumptions for shell thickness. Poisson’s ratio is the ratio of axial to lateral strain and can be extracted easily from a simple simulation of CNTs under tension, which yields values of 0.2, close to the value of 0.19 seen for graphene [58]. This leaves e0 and h as the only two parameters with unknown values. This section provides a brief review of previous methods used to calculate shell thickness 47 and nonlocal constant values and then explains an optimization approach to finding suitable values for the above parameters. 3.3.1 Review of previous work Regarding values of shell thickness, some previous researchers have used the interlayer spacing of graphene sheets in a graphitic structure, i.e., 3.4 Angstroms as the CNT shell thickness [13, 59]. Another way of addressing the wall thickness issue has been to use values of bending rigidity and in plane stiffness obtained from MD or ab initio simulations together with their classical continuum definitions to calculate values of shell thickness for use in equivalent continuum models [58, 60]. The continuum definition of bending rigidity D for a cylindrical shell is: 3 212(1 ) EhD υ= − , (3.9) and the tension rigidity C is defined as: C Eh= . (3.10) The bending and tension rigidities can be directly extracted from measurable characteristics of CNTs (such as the bending and tension energies) using MD simulations of CNT bending and CNT tension, respectively. Once the values of C and D have been determined the ratio of bending rigidity to tension rigidity is used to find the CNT shell thickness: 212(1 ) Dh C υ= − . (3.11) The aforementioned approach has been used by Kudin et al. where a value of 0.89 is reported for the CNT shell thickness [60]. Once the shell thickness is calculated, 48 knowing the values of in-plane stiffness, bending rigidity and torsional stiffness, other properties such as the Young and shear moduli can be determined. Previous research on finding the correct values of the nonlocal constant e0 for CNTs is very limited and there still is no consensus on the value of e0 that should be used for CNTs. Most previous publications look at nonlocal effects for a range of different e0 values, and present results similar to Figure (2.3) [36, 61-65]. Some publications have resorted to using the value of e0=0.39 reported by Eringen [31, 66]. However, it is important to note that e0 is a material parameter, i.e., the adopted value of coefficient e0 depends on the crystal structure and the physics of the problem under investigation, and it should be determined for each material independently [41]. Some recent publications have proposed e0 values by comparing nonlocal models with molecular dynamics simulation results. Zhang et al. [34] used a nonlocal Timoshenko beam model to compare axial buckling results with the molecular mechanics simulations of Sears and Brata [67] and proposed a value of 0.82 for the nonlocal constant. Hu et al. [38] compared nonlocal shell models for dispersion of propagating waves in CNTs with results from molecular dynamics simulations and reported e0 values of 0.6 for transverse waves and 0.2-0.23 for torsional waves. Based on the large scattering of the different e0 values used in previous publications, further investigation of this important material property for CNTs seems necessary, and could result in better understanding of nonlocal effects in CNTs. 3.3.2 Calculation of e0 and h using optimization In the work presented here it is proposed to determine the values of shell thickness and nonlocal constant based on a non-linear least-square fitting by minimizing the Euclidean norm of the difference between the buckling torques obtained directly from the MD simulations and the nonlocal elasticity shell model. This approach will also be used to determine the shell thickness for the classical shell models. In order to reduce the number of parameters to be determined from the optimization analysis, the surface shear modulus (G.h) is obtained directly from the MD torsional simulations and a Poisson’s ratio of 0.2 is used based on the existing literature [56]. Therefore, the shell thickness h is considered as the only optimization variable for 49 the classical continuum shell models (equations (2.18) and (2.31)) and both h and e0 are considered as optimization variables for the nonlocal elasticity based models (equations (2.17) and (2.30)). Values of h and e0 obtained from the optimization analysis are shown in Table 3.2. The analysis yields a wall thickness of 0.81 and 0.75 angstroms based on the classical Timoshenko shell model for the zigzag and armchair CNTs, respectively. The nonlocal Timoshenko shell model predicts a wall thickness of 0.85 and 0.86 angstrom and a nonlocal constant of 0.85 and 0.61 for armchair and zigzag nanotubes, respectively. We note that the values of shell thickness calculated here using the nonlocal shell models are very close to the value of 0.89 angstroms previously reported by Kudin et al [60] (only a 5% difference is observed). The e0 values obtained here are similar to those previously reported [34]. The difference between the values of e0 obtained for zigzag and armchair nanotubes may be attributed to the chirality-dependent anisotropy of CNT structures. Table 3.2- Values of shell thickness (h) and nonlocal constant (e0) obtained from non- linear least square fitting of MD simulation results with classical shell and nonlocal shell models. Armchair h(Ǻ) e0 Residual Norm (nN2.nm2) Classical Timoshenko 0.75 4.11 Nonlocal Timoshenko 0.85 0.85 0.09 Nonlocal Donnell 0.85 0.79 0.1 Zigzag h(Ǻ) e0 Residual Norm (nN2.nm2) Classical Timoshenko 0.81 5.9 Nonlocal Timoshenko 0.86 0.61 0.04 Nonlocal Donnell 0.86 0.57 0.05 50 3.3.3 Comparison of classical and nonlocal models with MD simulation results Although the values of shell thickness obtained from classical models are close to those previously reported based on theoretical concepts [68], for CNTs with small radii (e.g. less than 1 nm) the classical shell model does not show the correct trend of the critical buckling torque with changing diameter (Figure 3.7). This fact is verified by the high residual norm in Table 3.2. Figure 3.7- Comparison of buckling torques from classical shell models with MD results for (a) zigzag and (b) armchair CNTs based on properties given in Table 3.2. (a) (b) (10,0) (12,0) (16,0) (20,0) (14,0) (5,5) (6,6) (12,12) (10,10) (8,8) 51 Compared to MD simulations, the classical model shows a lower rate of increase in buckling torque with increase in CNT diameter. That is why in Figure 3.7 the trend line corresponding to the MD simulation has a steeper slope. This leads to the following important observation. The optimization algorithm is trying to minimize the residual norm and this goal is achieved by choosing values of wall thickness h to fit the classical model to the MD simulation results of the mid-point in the data set. For the data set used in this optimization, the mid-points correspond to (14,0) for the zigzag CNTs (Figure 7(a)) and (8,8) for armchair CNTs (Figure 7(b)). So if a different set of data (for example the first or last three points in each graph) is used for the optimization a different value for the wall thickness is obtained. It can be concluded that the classical models are not only unable to show the correct trend of change in buckling torque with change in diameter, but are also unable to determine a global wall thickness value that could be used for all zigzag or armchair nanotubes regardless of their diameter size. In comparison, the nonlocal elasticity shell models provide a much better fit with MD results (Figures 3.8 and 3.9) over a range of diameters because of the very low residual norms obtained from the analysis as shown in Table 3.2. For CNTs with small radii (e.g. less than 1 nm) the nonlocal shell model shows the correct trend of the critical buckling torque with changing diameter (Figures 3.8 and 3.9). Unlike the classical model, all the data points in the nonlocal model can be very closely fit to their corresponding values obtained from MD simulations. Thus, in the nonlocal model, regardless of which data points are used in the optimization algorithm, the calculated value for the wall- thickness are always very similar to those obtained in Table 3.2. In fact, when the optimization algorithm is performed on the nonlocal model using different combinations of 4 out of the 5 points in each dataset, the resulting values for wall thickness and nonlocal constant are similar to those reported in Table 3.2, and the resulting nonlocal model predicts the fifth point in the dataset with great accuracy. This is one of the benefits of using the nonlocal model, i.e., a global thickness can be used for CNTs of all diameters. 52 (a) (b) Figure 3.8- Comparison of buckling torques from classical and nonlocal Timoshenko shell models with MD results for (a) zigzag CNTs and (b) armchair CNTs based on properties given in Table 3.2. 53 (a) (b) Figure 3.9- Comparison of buckling torques from classical and nonlocal Donnell shell model with MD results for (a) zigzag CNTs and (b) armchair CNTs based on properties given in Table 3.2. 54 An important observation is that if the global thickness calculated using the nonlocal model is used in the classical model, the classical model yields buckling torques that are higher than those of the nonlocal models. The difference between predicted buckling torques from classical and nonlocal models is significant for CNTs with small diameters (as much as 40%) but as the CNT diameters increase, the classical buckling torques tend to converge to those of nonlocal models, and MD simulation results and the differences become negligible. This trend was predicted in chapter 2 (Figure 2.3) where the ratio of classical to nonlocal buckling torques were plotted against the CNT diameters and it was seen that this ratio converges to one as the CNT diameters increase. Another important observation is that changing the aspect ratio of CNTs while keeping within the limits of shell-shape buckling does not affect their torsional properties (rows 3 to 5 of Table 3.1-b). In other words, if the aspect ratio of the CNT changes within the limits of equation (3.8), the buckling torque, surface shear modulus and the buckling shear strain remain unchanged and the same set of values can be used for the nonlocal constant and shell thickness. Thus, the value of the nonlocal constant is independent of the magnitude of the geometric variables of the system. 55 3.4 MD simulations of Axial Buckling A question that remains is whether the values of nonlocal constant and shell thickness derived for the torsional buckling of CNTs (see section 3.3) can be applied to other loading conditions. This section attempts to answer this question by studying the axial buckling of CNTs. Previously, Zhang et al. [34] used a nonlocal Timoshenko beam model to compare axial buckling results with the molecular mechanics simulations of Sears and Brata [67] and proposed a value of 0.82 for the nonlocal constant, which shows good compatibility with the nonlocal constant values we have derived for torsional buckling of armchair CNTs. However Zhang’s results are obtained by comparing nonlocal models for axial buckling of CNTs with MD simulations of the beam-type axial buckling of CNTs. It would be of interest to study the shell-type axial buckling and see if the same values can be used. In this work shell-type axial buckling of several armchair CNTs is simulated using molecular dynamics and compared with the classical and nonlocal models seen in equations (2.40) and (2.39), respectively, in order to validate the wall thickness and nonlocal constant values obtained in the previous sections. 3.4.1 Axial buckling simulation setup Once again, quasi-static molecular dynamics is used to simulate the shell-type axial buckling (n >1) of several armchair nanotubes. The axial strain is simulated by mathematically changing the coordinates of the carbon atoms to a compressively strained state (Figure 3.10). These new coordinates are then input into the molecular dynamics simulator and the positions of the atoms at both ends of the CNT are fixed to simulate pined-pined boundary conditions. The CNT is then allowed to relax to an equilibrium configuration under the predefined axial strain. The number of relaxation steps required to reach equilibrium depends on the number of atoms in the system, the type of nanotube, the timestep and the applied strain. For each carbon nanotube the above simulation is performed for different compressive axial strains, and it is seen that above a certain value of strain, which is identified as the critical strain for axial buckling, the nanotube collapses into a buckled mode-shape when allowed to relax for a sufficient amount of time (Figure 3.11). 56 Figure 3.10- A compressive axial strain of 0.04 is induced in a (10,10) armchair CNT by mathematically changing the coordinates of carbon atoms (atoms at both ends of CNT shown in dark color are fixed). Figure 3.11- Different instances of relaxation of a (10,10) armchair nanotube under an axial compressive strain of 0.04. (a) simulation start (b) after relaxation for 19250 timesteps (c) after relaxation for 20250 timesteps (b) after relaxation for 30000 timesteps Axial compressive strain is induced 10 nm 9.6 nm 57 3.4.2 Axial buckling simulation results First, the axial buckling mode-shape calculated with molecular dynamics is checked against the mathematical assumption for the buckling mode-shape (equation (2.41)) to see if it is correct. As an example, based on equation (2.41), a (10,10) armchair nanotube with radius a=6.7 angstroms and length l=163.7 angstroms should have a buckling mode-shape with circumferential wave number n=2 and a longitudinal half- wave number m=5.2 corresponding to a longitudinal wave number of m/2=2.6. The buckled mode shape of a (10,10) armchair nanotube at critical strain of 0.035 simulated using molecular dynamics is shown in Figure 3.12. It is evident that there are two circumferential wavelengths in the transverse direction so that n=2. The number of longitudinal waves can also be counted with a visual inspection of Figure 3.12 and is close to m/2=2.6 as predicted by the analytical solution. Figure 3.12- Axial buckling mode-shape of (10,10) armchair nanotube at critical strain Once gain the progression of the potential energy of the system with strain (Figure 3.13) is used to find the compressional stiffness and values of surface Young’s modulus. These parameters are needed for use in equations (2.39) and (2.40). The surface Young’s modulus is easily found through the following relation: 58 2 2 1 2 d UEh al dπ ε= (3.12) Figure 3.13- Progression of potential energy with axial buckling strain of a (10,10) armchair nanotube. The critical axial buckling strains and the corresponding critical loads along with the surface Young’s moduli of several armchair nanotubes obtained from MD simulations are presented in Table 3.3. Table 3.3- Properties and critical forces of several armchair CNTs. Chiral Indices Diam.(Ǻ) Length (Ǻ) L / Diam E.h(GPa nm) Fcr (N.m) λ cr (6,6) 8.0 98.21 12.28 323.31 5.68E-08 0.84 (8,8) 10.7 131.00 12.24 307.93 5.41E-08 0.72 (10,10) 13.4 163.66 12.21 301.90 5.31E-08 0.65 (12,12) 16.1 165.00 10.25 299.56 5.26E-08 0.59 (14,14) 18.8 142.86 7.58 298.03 5.24E-08 0.54 (16,16) 21.5 160.00 7.44 296.10 5.20E-08 0.51 (20,20) 26.9 203.84 7.57 295.84 5.20E-08 0.46 Buckling 59 3.4.3 Values of nonlocal constant and shell thickness Once again, an optimization is performed to calculate the values of parameters to be used in the nonlocal and classical models. The shell thickness h is the optimization variable for the classical case and both shell thickness h and nonlocal constant e0 are optimization variables for the nonlocal model. The values obtained for these parameters are shown in Table 3.4. Table 3.4- Values of shell thickness (h) and nonlocal constant (e0) obtained from non- linear least square fitting of MD simulation results for axial buckling h(Ǻ) e0 Residual Norm (nN2) Classical Timoshenko 0.66 196.27 Nonlocal Timoshenko 0.81 0.94 6.66 Looking at the nonlocal model, the optimized values for shell thickness derived from the axial buckling simulations are consistent with those derived for torsional buckling in the previous section (less than 5% difference) while this is not true for the classical model (about 15% difference). The value of the nonlocal constant derived here is larger than that derived for torsional buckling and Zhang’s results for axial buckling [31] by about 10%; however, it is within the same range. This difference could be attributed to the different physics of the axial and torsional buckling problems, the differences in MD simulation setup (such as temperature control), as well as to small errors that exist in the numerical simulations. Once again it is seen that for smaller CNTs, the classical model is unable to show the correct rate of change in critical axial buckling loads with change in diameter (Figure 3.14 (a)) while the nonlocal shell model gives a much better prediction of molecular dynamics simulation results (Figure 3.14 (b)). It is seen that the axial buckling stresses from MD simulations are smaller than those predicted by the classical models. In comparison, the nonlocal model along with the thickness and nonlocal constants derived in the previous section give results that are a better match to MD simulations. The best match between MD simulations and nonlocal formulations is achieved for a nonlocal constant value of e0=0.94. The results from MD simulations, classical and nonlocal models are compared in Figure 3.14. 60 4 6 8 10 12 14 30 35 40 45 50 55 60 CNT Radius (Angstrom) A xi al B uc kl in g Lo ad (n N ) MD simulation Classical Model h=0.66 (a) 4 6 8 10 12 14 30 35 40 45 50 55 60 CNT Radius (Angstrom) A xi al B uc kl in g Lo ad (n N ) MD simulation Classical Model h=0.8 Nonlocal Model h=0.8 & e 0 =0.94 (b) Figure 3.14- Comparison of axial buckling loads from (a) classical and (b) nonlocal shell model with MD results for several armchair CNTs based on properties given in Table 3.3. e0 61 Chapter 4 NONLOCAL SHELL MODEL FOR TORSIONAL VIBRATION OF CNTS Torsional vibrations of CNTs are important to the design of different CNT based torsional devices including but not limited to sensors, oscillators and torsional pendulums. It is therefore necessary to have a clear understanding of CNT behavior under torsional vibrations. In this chapter a study of the dynamic torsional properties of CNTs is presented and dispersion relations of torsional waves are derived using a nonlocal shell model and compared with those of classical models and results from molecular dynamics simulations. The results presented are also compared with those derived in chapter 3 to validate the accuracy and efficiency of the nonlocal shell model compared to classical models. The results are also compared with those of Hu et al. [38], which, to the best of the author’s knowledge, are currently the most complete and up to date work regarding nonlocal models and values of nonlocal constant. Based on these comparisons a discussion is presented on the correct values of shell thickness and nonlocal constant and how the value of nonlocal constant may vary under different loading conditions. 4.1 Nonlocal Shell Model for CNT Torsional Waves and Vibrations Assuming a CNT can be modeled as a cylindrical shell with a uniform thickness (Figure 2.1) such that under torsion no in-plane elastic deformation or warping occurs, the displacements of any point in the cylindrical shell are [43]: 0 ; ; 0u v a wφ= = = , (4.1) where ( , )x tφ φ= is the rotation of the cross-section at x around the cylinder axis at time t, and this rotation depends only on the longitudinal coordinate x and is independent of θ. Thus the derivatives of all displacements with respect to θ are zero, i.e.: 62 ( ) 0 : , ,d u v w dθ Χ = Χ . (4.2) Adding inertial terms to equations (2.11) of chapter two and noting that u = w = 0, the following equilibrium relation holds for the torsional vibrations of a cylindrical shell: xdN Q q hv dx a θ θ θ ρ+ + =  , (4.3) where xdM Q dx θ θ= and qθ is an externally applied force. Thus [43]: x xdN dM q hv dx adx θ θ θ ρ+ + =  . (4.4) According to equations (2.9) and (2.10) of chapter 2: 2 (1 ) 2 2(1 )x R x K dv Eh dvN N dx dxθ θ υξ υ − ⎛ ⎞ ⎛ ⎞− ∇ = =⎜ ⎟ ⎜ ⎟+⎝ ⎠ ⎝ ⎠ , 3 2 (1 ) 2 12 (1 )x R x D dv EhM M a dx aθ θ υξ υ − ⎛ ⎞− ∇ = =⎜ ⎟ +⎝ ⎠ , (4.5) where 2/1K Eh υ= − is the membrane stiffness, 3 2/12(1 )D Eh υ= − is the bending rigidity and 20( )e dξ = is the nonlocal parameter. Substitution of equation (4.5) into (4.4) yields: 2 2 2 2 2(1 ) (1 )( )2(1 ) 12 R Eh h d v hv q a dx θ ξ ρυ + = − ∇ −+  . (4.6) 63 Now taking into account that 2 2/12 1h a  and introducing the shear wave speed /c G ρ= , in the absence of externally applied forces the above equation is simplified to: 2 2 2 2 2 0 d dc dx dx φ φφ ξ− − = . (4.7) This is the nonlocal form of the governing differential equation for the torsional vibrations of a cylindrical shell. In the absence of nonlocal effects 0ξ = and the above equation reduces to its classical form below [43]: 2 2 2 0 dc dx φφ − = . (4.8) 4.1.1 Propagating wave solution It can be assumed that a torsional wave of the form ( )0 i kx te ωφ φ −= is propagating along the axis of the cylindrical shell, where 2 /k π λ= is the wave number, λ is the wavelength and ω is the circular frequency of the propagating wave. For the plane wave mentioned above the following identities can be verified: 2φ ω φ= − , 2 2 2 d k dx φ φ= − , 2 2 2 2 d k dx φ ω φ= . (4.9) Substituting the above relations into equation (4.7) the following characteristic equation is obtained: 64 2 2 2 2 2( ) 0k c kω ξω φ− − + = . (4.10) Thus, for a torsional wave propagating in a cylindrical shell, the dispersion relation based on nonlocal elasticity is: 2 2 2 2 2 2 2 2 2 01 1 c k c k k e d k ω ξ= =+ + . (4.11) Once again, in the absence of nonlocal effects ( 0ξ = ), the above equation reduces to its classical form: 2 2 2c kω = . (4.12) It is interesting to look at and compare the group and phase velocities of propagating waves based on the nonlocal and classical solutions. I) Classical Solution The phase velocity Cp is the ratio of the circular frequency ω to the corresponding wave-number k. The group velocity Cg is the ratio of infinitesimal changes in ω and k (or the slope of the ω−k curve). The classical dispersion relation seen in equation (4.12) is used to obtain phase and group velocities: pC ck ω= = , g dC c dk ω= = . (4.13) It is seen that in the classical solution the phase and group velocities are equal to the shear speed of sound and are constant for all wave-numbers / frequencies, indicating no dispersion. 65 II) Nonlocal Solution Here the nonlocal dispersion relation seen in equation (4.11) is used to find group and phase velocities. The phase velocity is: 2 0.5 2 2 2 0.5 0(1 ) (1 ) p c cC k k e d k ω ξ= = =+ + , (4.14) and the group velocity is: 2 1.5 2 2 2 1.5 0(1 ) (1 ) g d c cC dk k e d k ω ξ= = =+ + . (4.15) It is seen that in the nonlocal solution, the group and phase velocities depend on the frequency/wave-number of the propagating wave, indicating dispersion. The strength of dispersion is governed by the positive nonlocal parameter ξ. The phase and group velocities are smaller in the nonlocal case than the classical value c and they decrease as the frequency of the waves increase. With increasing frequency, group velocity decreases at a faster rate compared to phase velocity. It is obvious to verify that classical results (equation (4.13)) ensue from the nonlocal results in the limiting case of 0ξ = . 4.1.2 Standing wave / vibration solution In the ensuing sections covering MD simulations of the dynamic torsional response of CNTs, it is seen that it is much easier to simulate and analyze the free vibrations (standing or non-propagating waves) of a CNT in order to see the real dispersion behavior. Therefore, the free torsional vibrations of a pined-pined CNT are simulated and the corresponding natural frequencies are studied to calculate group and phase velocities and analyze the dispersion behavior. Thus, in this section a brief discussion on the free torsional vibrations of a pined-pined CNT is presented. If ( , )x tφ is the rotation of the cross-section of the nanotube at a distance x from the support at time t the pined-pined boundary conditions state that: 66 (0, ) 0 , ( , ) 0t l tφ φ= = , (4.16) where l is the length of the nanotube. For each simulation, the nanotube is pre-deformed to a natural mode-shape at t=0 and released with zero initial velocity. Thus the initial conditions (at t=0) are: ( ,0) sin , 1,2,... , ( ,0) 0k x dx A k x l dt πφ φ′ ′= = = , (4.17) where k′ is the number of half waves along the length of the nanotube and l is the length of the nanotube. If ω is the circular frequency of the vibration it can be shown that: ( , ) sin cosk xx t A t l πφ ω′= , (4.18) where equation (4.18) is a solution of the governing differential equation (4.7) and satisfies the boundary and initial conditions of (4.16) and (4.17). Substituting this solution into equation (4.7) the following relation is obtained: 2 2 2 2 2 2 2( ) 0l k k cξ π ω π φ⎡ ⎤′ ′+ − =⎣ ⎦ . (4.19) The dispersion relation is found by equating the term in the brackets to zero: 2 2 2 2 2 2 2 22 2 2 2 2 2 2 02 2 ( ) ( ) 1 ( ) 1 ( ) k kc c l l k ke d l l π π ω π πξ ′ ′ = =′ ′+ + . (4.20) For comparison the classical dispersion relation is obtained by setting the nonlocal parameter to zero: 67 2 2 2 2 2( ) kc l πω ′= . (4.21) Taking note that k′ is related to the wave number k (k=k′π/l), it is seen that as expected, the dispersion relations of equation (4.20) and (4.21) are the same as the dispersion relations derived in equation (4.11) and (4.12), respectively. So the classical and nonlocal group and phase velocities derived in equations (4.13) to (4.15) hold for the case of free vibrations as well. This was anticipated from the beginning since free vibrations can be thought of as the superposition of two propagating waves of the same velocity and amplitude traveling in opposite directions. 4.1.3 Comparison of classical and nonlocal models Based on equations (4.12) and (4.21) classical models predict that for torsional waves (or torsional vibrations) in a cylindrical shell there is a linear relationship between the wave number (or modeshape) and the corresponding natural frequency. Classical models also predict that the phase and group velocities of torsional waves in a cylindrical shell are constant, independent of the wavenumber and equal to the shear speed of sound. In contrast to classical models, the nonlocal models derived in equations (4.11) and (4.20) show a nonlinear relation between the wave number of a torsional wave and its natural frequency. For higher wave numbers with smaller wave lengths, nonlocal models predict natural frequencies that are smaller than those predicted by classical models. This increased compliance of nonlocal models is due to the additional deformation mechanisms allowed in accounting for the total strain energy. A useful analogy to recall is the higher compliance of a Timoshenko beam model compared to Euler-Bernoulli beam model. It is well known that the Timoshenko model predicts lower values for natural frequencies and is indeed unavoidable for describing vibrations of thick beams where shear deformation is important, while the Euler- Bernoulli beam theory does not account for shear deformation. Nonlocal models predict that torsional waves in cylindrical shells are dispersive; as the wave numbers increase the corresponding phase and group velocities decrease. The highest phase and group velocities obtained from the nonlocal model correspond to the 68 first mode shape or the smallest wave number and are equal to the shear speed of sound. It is important to note that based on the nonlocal model, the magnitude of dispersion depends on the value of the nonlocal parameter ξ. It is interesting to investigate the general behavior of the nonlocal model, compare its nonlinearity to the linear classical model and quantify the relative errors that occur if classical models are used instead of the nonlocal models. Looking at the nonlocal relation between ω and k seen in equation (4.11) the Taylor series expansion can be written in the following form: 2 2 2 2 1/2 ( ) ( )2 3 2 5 ( ) (1 ) 1 1/ 2 3 / 4 . 1! 2! k k k c k ck k k ck c k c k ω ω ξξ ω ξ ξ −= → = ++ −⇒ = + + +" (4.22) It is seen that a polynomial describing the ω-k behavior should only have odd powers of k. The above expression can be used to find the general behavior of group and phase velocities: ( ) 2 2 4 ( ) 2 2 4 3 / 2 15 / 4 1! 2! 1/ 2 3 / 4 . 1! 2! k g k p d C c c k c k dk C c c k c k k ω ξ ξ ω ξ ξ −= = + + + −= = + + + " " (4.23) It is seen that polynomials describing the Cp-k and Cg-k behavior should only have even powers of k and are quadratic in nature. The relative errors in calculating ω(k), Cg and Cp using the classical model are: 69 ( ) ( ), ( ),( ) ( ), 2 2 4 2 2 4 1/ 2 3 / 4[1 ] 1! 2! 1/ 2 3 / 4 1! 2! k classical k nonlocal k k classical error ck ck k k ck k k ω ωω ω ξ ξ ξ ξ −= −− + + + = = − + " " ( ) , , , 2 2 4 2 2 4 3 / 2 15 / 4 1! 2! 3 / 2 15 / 4 1! 2! g classical g nonlocal g g classical C C error C C c c c k c k c k k ξ ξ ξ ξ −= −⎛ ⎞− + + +⎜ ⎟⎝ ⎠= = − + " " ( ) , , , 2 2 4 2 2 4 1/ 2 3 / 4 1! 2! 1/ 2 3 / 4 1! 2! p classical p nonlocal p g classical C C error C C c c c k c k c k k ξ ξ ξ ξ −= −⎛ ⎞− + + +⎜ ⎟⎝ ⎠= = − + " " (4.24) It is seen that the relative error is larger at higher wavenumbers. As an example, for a modeshape with a wavenumber of 1.3e10 corresponding to a wavelength of 4.8 angstroms, using a nonlocal constant value of e0=0.2 and an interatomic spacing of d=1.41 angstroms (ξ=e02d2), based on equation (4.24) the relative errors in calculating ω(k), Cg and Cp using the classical model are: ( ) ( ) ( ) 0.08 0.009 0.071 0.24 0.045 0.195 0.08 0.009 0.071 g p error error C error C ω = − + = − + = − + " " " For larger values of the nonlocal constant, the errors could be even more significant. The above comparisons can be better understood using a graphical 70 demonstration. Figure 4.1 (a) compares classical and nonlocal models by looking at the ratio of phase velocity to the shear speed of sound for various values of the wave number and different values of the nonlocal constant e0. 0 0.5 1 1.5 2 x 1010 0 0.2 0.4 0.6 0.8 1 1.2 k (1/m) C p/ c Classical Nonlocal e 0 =0.2 Nonlocal e 0 =0.4 Nonlocal e 0 =0.6 Nonlocal e 0 =0.8 (a) 0 0.5 1 1.5 2 x 1010 0 0.2 0.4 0.6 0.8 1 1.2 k (1/m) C g/ c Classical Nonlocal e 0 =0.2 Nonlocal e 0 =0.4 Nonlocal e 0 =0.6 (b) Figure 4.1- Comparison of classical and nonlocal models at different wavenumbers and for different values of the nonlocal constant e0 (a) phase velocities (b) group velocities 71 Figure 4.1 (b) shows the same comparison but for group velocities. The dispersive behavior of torsional waves in the nonlocal model and the non-dispersive behavior of torsional waves in the classical model are clearly observed in Figure 4.1. 4.2 MD simulations of Torsional Vibrations In order to observe the real dispersion behavior of torsional waves in CNTs and portray the superiority of nonlocal models compared to classical models, molecular dynamics is used to investigate the torsional vibration and torsional wave propagation characteristics of two armchair nanotubes. The purpose is to find values for the nonlocal constant and see how they match against those previously obtained in the third chapter. A comparison with the torsional wave propagation properties reported by Hu et al. [38] is also presented. The approach used previously [38] has been to study the propagation of torsional waves of different wavelengths and to calculate the corresponding phase velocities by measuring the time it takes a wave of certain wave length / frequency to travel between two designated points in a CNT. This approach is useful to some extent, but is limited in the sense that only a limited number of wavelengths can be studied, and in a MD simulation of a CNT, inducing torsional wave packets that purely consist of one certain frequency is a challenge by itself. In the present work it is proposed to simulate the free vibrations of pined-pined CNTs using molecular dynamics. Studying the free vibrations has several benefits. For example, it is possible to study a large number of the natural frequencies / mode-shapes and obtain a data set with sufficient data that can be used to draw confident conclusions. Also, once the natural frequency and mode-shape data is available, quantities such as the phase and group velocities can easily be calculated, helping answer important questions regarding the wave propagation and dispersion behavior of torsional waves in CNTs. 4.2.1 MD simulation setup In the present work molecular dynamics is used to simulate free torsional vibrations on (6,6) and (10,10) armchair nanotubes. Different natural mode-shapes of each nanotube are simulated and the corresponding natural frequencies are extracted 72 using the simulation outputs. To find the natural frequency corresponding to each mode- shape the nanotube is pre-deformed to that mode-shape by mathematically changing the coordinates (Figure 4.2) and rotating each row of atoms around the CNT axis according to the following relation: ( ,0) sin k xx A l πφ ′= , (4.25) where k′ is the mode-number (or half wave-length) corresponding to a certain mode- shape and l is the length of the CNT. Then the atoms at both ends are fixed to simulate pined-pined boundary conditions (Figure 4.2) and the CNT is allowed to vibrate freely for 10000 time steps of 0.5 fs. Mathematical transformation of atomic coordinates to the eight mode-shape Pinned end Pinned end Figure 4.2- Free torsional vibrations of a (6,6) armchair CNT at the eighth natural modes- shape corresponding to k′=8, where the mode-shape is induced by mathematically changing the atomic coordinates The rotation of the cross-sections of the CNT is extracted as an output of the simulation and plotted against time using MATLAB software. The natural vibration period is extracted by visualizing the free vibrations in the NanoHive post processor called Hive-keeper. The time period is then used to find the natural frequency of vibration corresponding to that mode-shape. Figure 4.3 shows the free torsional vibrations of a (6,6) armchair CNT at the eighth natural modes-shape corresponding to k′=8. At t=0 the 73 simulation starts and the nanotube is allowed to undergo free torsional vibrations. The rotation of the cross-sections of the CNT at different instances is shown. It is evident that the 8th mode-shape has a vibration period of τ =320 fs. t =0 t =27 fs t =54 fs t =80 fs t =106 fs t =133 fs t =160 fs Figure 4.3- Free torsional vibrations of a (6,6) armchair CNT at the 8th natural modes- shape corresponding to k′=8. Different instances of half a vibration period (τ /2=160 fs) are shown. An important consideration when performing the above molecular dynamics simulations is that the thermostat behaves as an active damper. By scaling the velocities to keep the temperature constant, the thermostat is randomly introducing or taking away kinetic energy to or from the system. Furthermore, the thermostatic kinetic energy is divided between atoms randomly and might excite mode-shapes other than the original mode-shape of the simulation. The purpose of the above simulations is to find the natural frequencies of pure mode-shapes in free torsional vibrations, and thus the thermostat which acts as an external stimulus is not applied to the system. Omitting thermostatic effects, the only change in kinetic energy is due to the conversion of potential energy of vibration and vise versa, and the total energy of the system remains constant. 74 The lengths of the CNTs were chosen such that each CNT has 145 rows of carbon atoms in the axial direction. The number of rows chosen above has the benefit that most of the simulated mode-shapes have nodes that coincide with the spatial position (longitudinal coordinate) of a certain row of atoms. Since the axial distance between two consecutive rows in an armchair nanotube is 1.21 angstroms, the total length of each simulated nanotube is 174.24 angstroms. 4.2.2 MD simulation results The torsional vibration simulations performed using molecular dynamics give the following results for the vibration periods and natural frequencies of the different mode- shapes of (6,6) and (10,10) armchair nanotubes, shown in Tables 4.1 and 4.2, respectively. Table 4.1- Vibration periods and natural frequencies of several different mode-shapes of (6,6) armchair CNT. Mode-shape k' Wavelength λ (m) Period τ (s) Frequency f (Hz) Angular Freq. ω (rad/s) 1 3.50E-08 2.52E-12 3.97E+11 2.50E+12 2 1.75E-08 1.26E-12 7.92E+11 4.97E+12 4 8.75E-09 6.36E-13 1.57E+12 9.88E+12 8 4.38E-09 3.20E-13 3.13E+12 1.96E+13 9 3.89E-09 2.86E-13 3.50E+12 2.20E+13 12 2.92E-09 2.10E-13 4.76E+12 2.99E+13 16 2.19E-09 1.60E-13 6.27E+12 3.94E+13 18 1.94E-09 1.41E-13 7.09E+12 4.46E+13 24 1.46E-09 1.07E-13 9.39E+12 5.90E+13 36 9.72E-10 7.10E-14 1.41E+13 8.85E+13 48 7.29E-10 5.40E-14 1.85E+13 1.16E+14 54 6.48E-10 4.80E-14 2.08E+13 1.31E+14 66 5.30E-10 4.00E-14 2.50E+13 1.57E+14 72 4.86E-10 3.68E-14 2.72E+13 1.71E+14 75 Table 4.2- Vibration periods and natural frequencies of several different mode-shapes of (10,10) armchair CNT. Mode-shape k' Wavelength λ (m) Period τ (s) Frequency f (Hz) Radial Freq. ω (rad/s) 1 3.50E-08 2.61E-12 3.83E+11 2.41E+12 4 8.75E-09 6.53E-13 1.53E+12 9.62E+12 16 2.19E-09 1.64E-13 6.11E+12 3.84E+13 24 1.46E-09 1.09E-13 9.17E+12 5.76E+13 36 9.72E-10 7.37E-14 1.36E+13 8.53E+13 54 6.48E-10 5.00E-14 2.00E+13 1.26E+14 66 5.30E-10 4.15E-14 2.41E+13 1.51E+14 72 4.86E-10 3.80E-14 2.63E+13 1.65E+14 4.3 Comparison of MD Simulation Results with Classical and Nonlocal Models Both the classical and nonlocal models for the radial frequency of torsional vibrations, seen in equations (4.21) and (4.20), respectively, are dependent on the value of the shear speed of sound (c) in CNTs which is equal to: . . G G hc hρ ρ= = , (4.26) where G.h is the surface shear modulus previously obtained in chapter 3, ρ is the mass density of the CNT and ρ.h is the surface mass density found using the following relation: . 2 2 m m mh V ahl al ρ ρπ π= = ⇒ = . (4.27) Thus, knowing the number of carbon atoms in the CNT, the mass of each carbon atom and the radius a and length l of the CNT, the surface mass density ρ.h is easily calculated and used in equation (4.26) along with the value of the corresponding surface 76 shear modulus G.h to find the shear speed of sound. The values of the shear speed of sound calculated for (6,6) and (10,10) armchair nanotubes using the above approach are shown in Table 4.3. Table 4.3- Surface shear modulus, surface density and shear wave speed c for (6,6) and (10,10) armchair nanotubes CNT G.h (N/m) ρ.h (kg/m2) c (m/s) (6,6) 153.7 7.88E-7 1.39E4 (10,10) 145.3 7.88E-7 1.35E4 4.3.1 Comparison of radial frequencies Based on the classical model of equation (4.21) and the values of sound speed given in Table 4.3 the radial frequencies of the first mode-shape (k′=1) of torsional vibrations for (6,6) and (10,10) armchair CNTs are: 1 1 (6,6) 2.51 12 ( / ) (10,10) 2.42 12 ( / ) k k E rad s E rad s ω ω ′= ′= → = → = These values are almost identical to the radial frequencies obtained from molecular dynamics shown in Tables 4.1 and 4.2 (the error is less than 0.5%). If the same classical model is used to find the radial frequencies of the 72nd mode-shape (k′=72) of torsional vibrations for (6,6) and (10,10) armchair CNTs the following values are obtained. 72 72 (6,6) 1.80 14 ( / ) (10,10) 1.74 14 ( / ) k k E rad s E rad s ω ω ′= ′= → = → = The radial frequency values obtained from the classical model are higher than those of MD simulation results (shown in Tables 4.1 and 4.2) and have errors of more 77 than 5%. It is evident that at these higher mode-shapes with smaller wavelengths, the classical model loses its accuracy (Figure 4.4), as discussed in section 4.1.3. 0 20 40 60 80 0 10 20 30 40 50 60 70 80 Mode-shape Number k ′ N or m al iz ed F re qu en cy ( ω / ω 1 ) MD simulation Classical Model Nonlocal Model e 0 =0.185 (a) 0 20 40 60 80 0 10 20 30 40 50 60 70 80 Mode-shape Number k ′ N or m al iz ed F re qu en cy ( ω/ ω 1) MD simulation Classical Model Nonlocal Model e 0 =0.18 (b) Figure 4.4- Comparison of local and nonlocal models for the normalized radial frequency versus mode-shape number with MD simulation results (a) (6,6) CNT (b) (10,10) CNT. 78 In comparison, looking at the nonlocal model of equation (4.20) the term in the denominator grows with increasing values of wave number and could compensate for the size effects that become significant at smaller wave-lengths. A least square optimization algorithm is used to find the correct value of the nonlocal parameter ξ (or nonlocal constant e0) to be used in the nonlocal model. The optimized values for the nonlocal constant are found as 0.185 and 0.18 for the (6,6,) and (10,10) nanotubes, respectively. It is seen that using these nonlocal constant values the nonlocal model can predict the MD simulation results with great accuracy, both at small and large wavelengths (Figure 4.4). The above results obtained for the nonlocal constant values are comparable to those previously reported. Hu et al. [38] compared nonlocal shell models for dispersion of torsional waves in CNTs with results from molecular dynamics simulations and reported an e0 value of 0.2 for (10,10) armchair CNTs. 4.3.2 Comparison of group velocities Based on equation (4.13) the group velocity is the slope of the natural frequency versus wave-number (ω-k) curve. To find the slope, a numerical differentiation technique can be used. For example, a two point numerical differentiation would look like: 1 ( 1) k kslope at wavenumber k k k ω ω −−= − − (4.28) When using the above numerical differentiation technique, small errors in reading the MD simulation output can lead to large errors once the derivation is carried out. A better way of finding the slope would be to use curve-fitting to find a polynomial that portrays the general behavior of the ω-k curve and then use mathematical differentiation to find the derivative of the polynomial and calculate its value at different wavenumbers. As the nonlocal model gives a precise fit to MD simulations, a polynomial based on the nonlocal model could accurately model the MD ω-k behavior. Based on equation (4.22) and the discussion presented in section 4.1.3, different polynomials are fit to the ω-k curve and it is seen that a third order polynomial (ω(k)=ak+bk3) can model the ω-k behavior with sufficient accuracy. This third order 79 polynomial is then differentiated to find the group velocities at different wavenumbers. The MD group velocities for (6,6) and (10,10) armchair nanotubes are presented in Table 4.4 and compared with the group velocities obtained from the classical and nonlocal models of equations (4.21) and (4.20), respectively. Table 4.4- Comparison of group velocities of classical and nonlocal models with MD simulation results for (a) (6,6) and (b) (10,10) armchair Mode-shape k' Wavenumber k (1/m) Cg MD (m/s) Cg Classical (m/s) Cg Nonlocal (m/s) 1 1.795E+08 1.390E+04 1.393E+04 1.393E+04 2 3.590E+08 1.390E+04 1.393E+04 1.392E+04 4 7.181E+08 1.389E+04 1.393E+04 1.392E+04 8 1.436E+09 1.387E+04 1.393E+04 1.390E+04 9 1.616E+09 1.387E+04 1.393E+04 1.389E+04 12 2.154E+09 1.384E+04 1.393E+04 1.386E+04 16 2.872E+09 1.380E+04 1.393E+04 1.381E+04 18 3.231E+09 1.377E+04 1.393E+04 1.378E+04 24 4.308E+09 1.367E+04 1.393E+04 1.367E+04 36 6.463E+09 1.339E+04 1.393E+04 1.335E+04 48 8.617E+09 1.299E+04 1.393E+04 1.293E+04 54 9.694E+09 1.275E+04 1.393E+04 1.269E+04 66 1.185E+10 1.219E+04 1.393E+04 1.215E+04 72 1.293E+10 1.186E+04 1.393E+04 1.190E+04 (a) Mode-shape k' Wavenumber k (1/m) Cg MD (m/s) Cg Classical (m/s) Cg Nonlocal (m/s) 1 1.795E+08 1.350E+04 1.348E+04 1.348E+04 4 7.181E+08 1.349E+04 1.348E+04 1.347E+04 16 2.872E+09 1.340E+04 1.348E+04 1.336E+04 24 4.308E+09 1.327E+04 1.348E+04 1.321E+04 36 6.463E+09 1.299E+04 1.348E+04 1.289E+04 54 9.694E+09 1.235E+04 1.348E+04 1.222E+04 66 1.185E+10 1.179E+04 1.348E+04 1.167E+04 72 1.293E+10 1.146E+04 1.348E+04 1.137E+04 (b) 80 The classical model predicts a constant group velocity equal to the shear speed of sound c. The nonlocal model predicts group velocities that decrease in a nonlinear manner with the increase of the corresponding wave numbers. Once again, it is seen that the classical model is accurate for smaller wave numbers (larger wavelengths) but loses its accuracy at higher wave numbers. If the classical model is used for higher wavenumbers it can result in errors of more than 15% when compared to MD simulations (Figure 4.5). In contrast to the classical model, the nonlocal model is accurate at all wave numbers (wavelengths) and when used along with the nonlocal constants derived in section 4.3.1 it precisely conveys the dispersion behavior seen in MD simulations. Figure 4.5 is based on the results of Table 4.4 and compares the classical and nonlocal models for group velocity with the values obtained from MD simulations. It is interesting to see how well the nonlocal model can predict the phase velocities obtained from MD simulations. This is done by studying the “goodness of fit” in figure 4.5 through the chi-square test. The chi-square parameter is the sum of the differences between observed and expected outcomes, each squared and divided by the expectation: 2 2 ( )O E E χ −= ∑ . (4.29) In the current example, the observed quantities are the MD simulation results and the expected outcomes are the results from the nonlocal model. Figures 4.5 (a) and 4.5 (b) have chi-square parameter ( χ2) values of 1.49e-4 and 2.26e-4, respectively. The χ2 values are very small compared to the average value of the observed or expected outcome (which is between 0.84 and 1), and this indicates a very high goodness of fit. 81 0 2 4 6 8 10 12 14 x 109 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 Wavenumber k (1/m) N or m al iz ed G ro up V el oc ity ( C g / C g1 ) MD simulation Classical Model Nonlocal Model e 0 =0.185 (a) 0 2 4 6 8 10 12 14 x 109 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 Wavenumber k (1/m) N or m al iz ed G ro up V el oc ity ( C g / C g1 ) MD simulation Classical Model Nonlocal Model e 0 =0.18 (b) Figure 4.5- Comparison of local and nonlocal models for the normalized group velocity versus wavenumber with MD simulation results (a) (6,6) CNT (b) (10,10) CNT. 82 4.3.3 Comparison of phase velocities The phase velocities corresponding to MD simulations can easily be found by finding the ratios of the natural frequencies ω to their corresponding wave numbers k. The polynomial expression found in the previous section describing the ω-k behavior of MD simulations is used to find the ω/k ratio and calculate the phase velocity corresponding to each wavenumber k. The phase velocity values of the different modeshapes of (6,6) and (10,10) armchair CNTs obtained from molecular dynamics are shown in Table 4.5 along with the phase velocities calculated using the classical and nonlocal models (equations (4.13) and (4.14)). Table 4.5- (a) Comparison of phase velocities of classical and nonlocal models with MD simulation results for (6,6) armchair CNT Mode-shape k' Wavenumber k (1/m) Cp MD (m/s) Cp Classical (m/s) Cp Nonlocal (m/s) 1 1.795E+08 1.390E+04 1.397E+04 1.397E+04 2 3.590E+08 1.390E+04 1.397E+04 1.397E+04 4 7.181E+08 1.390E+04 1.397E+04 1.396E+04 8 1.436E+09 1.389E+04 1.397E+04 1.396E+04 9 1.616E+09 1.389E+04 1.397E+04 1.395E+04 12 2.154E+09 1.388E+04 1.397E+04 1.395E+04 16 2.872E+09 1.387E+04 1.397E+04 1.393E+04 18 3.231E+09 1.386E+04 1.397E+04 1.392E+04 24 4.308E+09 1.382E+04 1.397E+04 1.388E+04 36 6.463E+09 1.373E+04 1.397E+04 1.379E+04 48 8.617E+09 1.360E+04 1.397E+04 1.365E+04 54 9.694E+09 1.352E+04 1.397E+04 1.357E+04 66 1.185E+10 1.333E+04 1.397E+04 1.339E+04 72 1.293E+10 1.322E+04 1.397E+04 1.328E+04 83 Table 4.5- (b) Comparison of phase velocities of classical and nonlocal models with MD simulation results for (10,10) armchair CNT Mode-shape k' Wavenumber k (1/m) Cp MD (m/s) Cp Classical (m/s) Cp Nonlocal (m/s) 1 1.795E+08 1.350E+04 1.358E+04 1.358E+04 4 7.181E+08 1.350E+04 1.358E+04 1.358E+04 16 2.872E+09 1.347E+04 1.358E+04 1.354E+04 24 4.308E+09 1.342E+04 1.358E+04 1.349E+04 36 6.463E+09 1.333E+04 1.358E+04 1.338E+04 54 9.694E+09 1.312E+04 1.358E+04 1.314E+04 66 1.185E+10 1.293E+04 1.358E+04 1.294E+04 72 1.293E+10 1.282E+04 1.358E+04 1.283E+04 The classical model predicts a constant phase velocity equal to the shear speed of sound c. The nonlocal model predicts phase velocities that decrease in a nonlinear manner with the increase of the corresponding wave numbers. Once again it is seen that the classical model is accurate for smaller wave numbers (larger wavelengths) but loses its accuracy at higher wave numbers. In contrast to the classical model, the nonlocal model is accurate at all wave numbers (wavelengths) and used along with the nonlocal constants derived in section 4.3.1 it can precisely convey the dispersion behavior obtained from MD simulations. Figure 4.6 is based on the results of Table 4.5 and compares the classical and nonlocal models for phase velocity with the values obtained from MD simulations. It is important to note that the dispersion behavior of phase velocities seen in the nonlocal model as well as molecular dynamic simulations is similar to the behavior previously reported by Hu et al [38]. However, the phase velocity derived for the first mode-shape in the present work (which is equal to the shear speed of sound) is larger than the value previously reported by Hu et al. This is probably due to the higher torsional stiffness or surface shear moduli that has been derived for the armchair nanotubes in the present work, which in turn could be caused by the use of different potential functions in the MD simulations. The approach used for finding the dispersion relations in the present work is more accurate compared to that of Hu et al. due to several factors: 84 - Hu et al. have looked at the velocity of torsional waves propagating along a CNT using MD simulations of forced wave propagation. The method used has less control over the frequency range of the excited waves and it is difficult to make sure that an excited wave consists of a single phase/ single frequency. It is likely that waves excited in such manner would be wave packets with a range of frequencies. In comparison, the method of free vibrations used in the present work allows greater control over the excited frequencies. It is possible to visualize the modeshape over several hundred vibration periods and make sure that only one frequency is excited. - The number of data points reported by Hu et al. is limited, and most data points have been reported for the small frequency / wavenumber range. Based on the discussions presented in the current work, it is evident that nonlocal effects become significant at higher wavenumbers and so there is a need to have sufficient simulation data at higher frequencies. The approach presented in the current work allows the simulation of different frequencies at the high end of the frequency spectrum, providing a rich dataset of simulation results to compare with the classical and nonlocal models. - In the present work, the MD simulation outputs are vibration frequencies, which can be used to calculate both the group and the phase velocities of propagating waves, in contrast to the previous work, which only considers the phase velocities. - In the present work the AIREBO potential field has been used, which is currently one of the most accurate potential fields available for the MD simulation of CNTs. Hu et al. have used the second generation REBO potential. Regarding the accuracy of the nonlocal model in predicting the MD simulation results for phase velocities, Figures 4.5 (a) and 4.5 (b) have chi-square parameter ( χ2) values of 6.94e-6 and 7.457e-5, respectively, which once again indicate a very high degree of fit of the nonlocal model to the MD simulation results. 85 0 2 4 6 8 10 12 14 x 109 0.85 0.9 0.95 1 1.05 Wavenumber k (1/m) N or m al iz ed P ha se V el oc ity ( C p / C p1 ) MD simulation Classical Model Nonlocal Model e 0 =0.185 (a) 0 2 4 6 8 10 12 14 x 109 0.85 0.9 0.95 1 1.05 Wavenumber k (1/m) N or m al iz ed P ha se V el oc ity ( C p / C p1 ) MD simulation Classical Model Nonlocal Model e 0 =0.18 (b) Figure 4.6- Comparison of local and nonlocal models for the normalized phase velocity versus wavenumber with MD simulation results. (a) (6,6) CNT (b) (10,10) CNT 86 4.4 Discussion on the Nonlocal Constant Value As mentioned in chapter 1, one of the main goals of the current work is to shed light on the inconsistency of the values reported for the nonlocal constant e0 of CNTs. Previous literature which reports e0 values for CNTs is very limited. Zhang et al. [34] used a nonlocal Timoshenko beam model to compare axial buckling results with the molecular mechanics simulations of Sears and Brata [67] and proposed a value of 0.82 for the nonlocal constant. Hu et al. [38] compared nonlocal shell models for dispersion of propagating waves in CNTs with results from molecular dynamics simulations and reported e0 values of 0.6 for transverse waves and 0.2-0.23 for torsional waves. The study carried out in the present work shows e0 values of 0.6-0.8 for torsional buckling (chapter 3), e0 values of about 0.9 for axial buckling (chapter 3) and e0 values of about 0.18 for torsional vibrations (chapter 4). What is the reason behind this scattering of the obtained e0 values? All the e0 values mentioned above have been obtained by comparison of nonlocal models with results of molecular dynamics simulations. Previous molecular dynamics simulations have been performed using different potential functions which could lead to different results. Thus, the accuracy of the molecular dynamics simulations could affect the reported values for the nonlocal constant e0. It is important to emphasize once again that all the MD simulations performed in the present work were carried out using the AIREBO potential function. Another reason could be due to the different loading conditions and the physics involved in each one of the aforementioned problems. Eringen [41] states that the adopted value of coefficient e0 depends on the crystal structure and the physics of the problem under investigation, and it should be determined for each material independently [41]. Although all the reported values mentioned above are for one material with a specified atomic structure (single-walled carbon nanotube), they have been reported for different loading conditions, which alter the physics of each problem. Perhaps one of the most important factors which has not been previously reported is the directional dependency of the nonlocal constant. To better understand this matter the nonlocal formulas obtained for the torsional and axial buckling loads and the 87 frequency of torsional vibrations are compared to the corresponding classical models here: 2 2 2 2 0 21 ( ) classical nonlocal MTorsional Buckling M d m ae n la π→ = ⎛ ⎞+ +⎜ ⎟⎝ ⎠ 2 2 2 2 0 21 ( ) classical nonlocal FAxial Buckling F d m ae n la π→ = ⎛ ⎞+ +⎜ ⎟⎝ ⎠ 2 2 2 2 2 0 21 ( ) classical nonlocalTorsional Vibrations d k ae la ωω π→ = ′⎛ ⎞+ ⎜ ⎟⎝ ⎠ (4.30) It is seen that generally the relation between the classical and nonlocal values of a characteristic property P of a cylindrical shell modeling a CNT system is of the following form: 2 2 2 2 0 21 ( ) classical nonlocal PP d m ae n la π= ⎛ ⎞+ +⎜ ⎟⎝ ⎠ (4.31) It is seen that the nonlocal constant e0 which controls the magnitude of the nonlocal effect is multiplied by the terms in the bracket which are functions of the number of longitudinal waves (m or k′) and the number of circumferential waves (n). Based on which of the terms in the bracket is larger in magnitude, different values are obtained for e0. This is better understood by studying the e0 values in Table 4.6. 88 Table 4.6- Obtained values of nonlocal constant Problem definition Problem type Typical value of (mπa/l)2 Typical value of n2 Dominating factor Value of e0 Torsional buckling static 0.7 4 n2 0.6-0.8 Axial buckling static 0.6 4 n2 0.8-0.95 Torsional vibrations dynamic 5 0 (mπa/l)2 0.18-0.23 Transverse vibrations dynamic 0.2 1 n2 0.6 Based on the results given in Table 4.6, it is evident that whenever there are waves in the circumferential direction such that n2 > (mπa/l)2, n2 is the dominating factor controlling the magnitude of the term (mπa/l)2 + n2 and the nonlocal constant values obtained are in the 0.6-0.95 range. In contrast, whenever there are no waves in the circumferential direction or the number of waves in the circumferential direction is small compared to the number of waves in the longitudinal direction such that (mπa/l)2 > n2, (mπa/l)2 is the dominating factor controlling the magnitude of the term (mπa/l)2 + n2 and the nonlocal constant values obtained are in the 0.18-0.23 range. The above observation leads to the conclusion that e0 must be direction dependent. That is, based on the loading conditions and deformation mechanisms that take place in a specific problem, the existence of circumferential and/or longitudinal waves in the deformed shape of the CNT dictates the correct value of e0 that should be used in the nonlocal model. 89 Chapter 5 SUMMARY AND CONCLUSIONS 5.1 Summary of Present Work and Major Findings (1) A nonlocal elasticity model for the torsional buckling of single-walled carbon nanotubes (SWCNTs) is developed. Eringen’s nonlocal elasticity theory is used in conjunction with the classical Timoshenko and Donnell shell models to construct a modified nonlocal continuum model for torsional buckling that can be used for all SWCNTs regardless of their size, specifically their diameter. It is shown that classical models always predict the same value for the non-dimensional buckling load of CNTs, while nonlocal models predict non-dimensional buckling loads that depend on the geometric dimensions and aspect ratios of the CNT, allowing the incorporation of “small- size” effects. As the lengths l and radii a of CNTs become smaller and comparable to the inter-atomic distance d, the d/a and d/l ratios become larger and have a greater impact on the magnitude of buckling loads predicted by the nonlocal models. In addition to the geometric aspect ratios, the magnitude of buckling loads predicted by the nonlocal models also depends on the buckling mode-shape and the value of the nonlocal constant e0. The difference between the buckling loads predicted from classical and nonlocal models is significant for CNTs with small radii, but this difference becomes negligible at larger radii where the nonlocal and classical models converge. (2) The classical molecular dynamics method is used to simulate the torsional buckling of a variety of different CNTs, to extract mechanical properties such as the surface shear moduli and buckling torques, and to establish the validity of the nonlocal model over the classical models in predicting the buckling loads of CNTs. The values of shell thickness h and nonlocal elastic constant e0 are found by fitting the continuum models (classical and nonlocal) to MD simulation results based on the least-squares technique. The classical shell model yields a wall thickness of 0.81 and 0.75 angstroms and the nonlocal model yields a wall thickness of 0.85 and 0.86 angstrom and a nonlocal constant of 0.85 and 0.61 for armchair and zigzag nanotubes, respectively. For CNTs with small radii (e.g. less than 1 nm) the classical shell model does not show the correct 90 trend of the critical buckling torque with changing diameter and could overestimate the buckling torque by as much as 40%. In contrast, the nonlocal shell model predicts the correct value of the critical buckling torque at all diameters. An important observation is that the values of the nonlocal constant and shell thickness are independent of the magnitude of the geometric variables of the system. Selected MD simulations of axial buckling of SWCNTs are also presented. Using the same shell thickness values obtained for torsional buckling, the nonlocal model for axial buckling is fit to the MD simulation results and a slightly higher e0 value is obtained, which may be due to the different loading conditions and deformation mechanisms that take place. (3) The dispersion of torsional waves due to nanoscale size-effects is studied for pinned-pinned SWCNTs. It is seen that the classical model is non-dispersive and predicts constant phase and group velocities while the nonlocal model is dispersive and predicts phase and group velocities that depend on the wavelength and the value of the nonlocal constant e0. The MD results for the free vibrations of pinned-pinned (6,6) and (10,10) armchair CNTs are fitted to the nonlocal model to extract the value of the nonlocal constant which is found to be approximately 0.18. Using this value of the nonlocal constant, the nonlocal model can predict the phase and group velocities of torsional waves of all wavelengths with great accuracy, while the classical model loses its accuracy at smaller wavelengths. If the classical model is used for torsional waves with wavelengths smaller than 8 angstroms (which is about 5 times the interatomic spacing d) it could overestimate the group and phase velocities by as much as 15% and 5%, respectively, and this error is significant for wavelengths of up to 20 angstroms. (4) Results obtained for torsional buckling, axial buckling, torsional vibrations and transverse vibrations are compared and it is proposed that the inconsistency of the values reported for the nonlocal constant in previous literature is partially due to directional dependency. That is, whenever there are circumferential waves in the deformed shape of a CNT, the circumferential wavenumber dominates the behavior of the nonlocal model and the values obtained for the nonlocal constant are in the 0.6 to 0.95 range. In the absence of circumferential waves in the deformed shape of a CNT, the longitudinal wavenumber dominates the behavior of the nonlocal model and the values obtained for the nonlocal constant are between 0.18 and 0.23. 91 5.2 Suggestions for Future Work Based on the importance of nonlocal continuum modeling to the accurate design and characterization of nanoscale devices and the findings of this thesis, it is recommended that the following studies be undertaken; (1) In the present work, torsional and axial buckling loads are derived using the equilibrium of forces approach. The equilibrium approach is useful for studying properties at the critical buckling strains/loads where the buckling deformations are assumed to be small; however, this method cannot be used to study the post buckling behavior where buckling deformations are quite large. Based on the results obtained in the present work, it is seen that size-effects are very significant at higher buckling modeshapes and so the post buckling region could be the subject of a further study. It is therefore recommended that the energy approach be used to develop a modified nonlocal elasticity based model and study size-effects in the post buckling region. (2) Based on the limitations of the molecular dynamics simulation software used in the present work, all simulations are based on CNTs with pinned-pinned boundary conditions, and so the nonlocal models presented have been developed for pinned-pinned boundary conditions . Modification of the molecular dynamics simulator and addition of features that allow the consideration of other boundary conditions could yield a rich set of data, useful for the extension and further validation of the nonlocal models presented in the current work. 92 Bibliography [1] Taniguchi, N., 1974, "On the Basic Concept of Nanotechnology," Proceedings of the International Conference of Production Engineering, pp. 18–23. [2] Drexler, K.E., 1992, Nanosystems: molecular machinery, manufacturing, and computation, John Wiley & Sons Inc., New York, NY, USA, pp. 1-576. [3] Capek, I., 2006, Nanocomposite structures and dispersions: science and nanotechnology-fundamental principles and colloidal particles, Elsevier, New York, NY, USA, pp. 1-293. [4] Iijima, S., 1991, "Helical Microtubules of Graphitic Carbon," Nature, 354(6348), pp. 56-58. [5] Zheng, L. X., O'Connell, M. J., Doorn, S. K., Liao, X. Z., Zhao, Y. H., Akhadov, E. A., Hoffbauer, M. A., Roop, B. J., Jia, Q. X., and Dye, R. C., 2004, "Ultralong Single-Wall Carbon Nanotubes," Nature Materials, 3(10), pp. 673-676. [6] Reich, S., Thomsen, C., and Maultzsch, J., 2004, Carbon nanotubes: basic concepts and physical properties, Wiley-VCH, pp. 1-224. [7] Kawasaki, S., Kanamori, Y., Iwai, Y., Okino, F., Touhara, H., Muramatsu, H., Hayashi, T., Kim, Y. A., and Endo, M., 2008, "Structural Properties of Pristine and Fluorinated Double-Walled Carbon Nanotubes Under High Pressure," Journal of Physics and Chemistry of Solids, 69(5-6), pp. 1203-1205. [8] Tans, S. J., Devoret, M. H., Dai, H., Thess, A., Smalley, R. E., Geerligs, L. J., and Dekker, C., 1997, "Individual Single-Wall Carbon Nanotubes as Quantum Wires," Nature, 386(6624), pp. 474-477. [9] Martel, R., Schmidt, T., Shea, H. R., Hertel, T., and Avouris, P., 1998, "Single-and Multi-Wall Carbon Nanotube Field-Effect Transistors," Applied Physics Letters, 73(17), pp. 2447-2449. [10] Postma, H. W. C., Teepen, T., Yao, Z., Grifoni, M., and Dekker, C., 2001, "Carbon Nanotube Single-Electron Transistors at Room Temperature," Science, 293(5527), pp. 76-79. [11] Gotou, Y., Matsumoto, K., Yasutake, M., Muramatsu, H., and Kim, J. M., 2002, "Single Wall Carbon Nanotube Cantilever: Fabrication and Application," Journal of the Surface Science Society of Japan, 23(2), pp. 116-122. [12] Jensen, K., Kim, K., and Zettl, A., 2008, "An Atomic-Resolution Nanomechanical Mass Sensor," Nature Nanotechnology, 3(9), pp. 533-537. 93 [13] Hall, A. R., An, L., Liu, J., Vicci, L., Falvo, M. R., Superfine, R., and Washburn, S., 2006, "Experimental Measurement of Single-Wall Carbon Nanotube Torsional Properties," Physical Review Letters, 96(25), pp. 256102. [14] Hall, A. R., Falvo, M. R., Superfine, R., and Washburn, S., 2007, "Electromechanical Response of Single-Walled Carbon Nanotubes to Torsional Strain in a Self-Contained Device," Nature, 2, pp. 413-416. [15] Hall, A. R., Falvo, M. R., Superfine, R., and Washburn, S., 2008, "A Self-Sensing Nanomechanical Resonator Built on a Single-Walled Carbon Nanotube," Nano Letters, 8(11), pp. 3746-3749. [16] Bourlon, B., Glattli, D. C., miko, c., Forro, L., and Bachtold, A., 2004, "Carbon Nanotube Based Bearing for Rotational Motions," Nano Letters, 4(4), pp. 709-712. [17] Papadakis, S. J., Hall, A. R., Williams, P. A., Vicci, L., Falvo, M. R., Superfine, R., and Washburn, S., 2004, "Resonant Oscillators with Carbon-Nanotube Torsion Springs," Physical Review Letters, 93, pp. 146101-146101. [18] Cohen-Karni, T., Segev, L., Srur-Lavi, O., Cohen, S. R., and Joselevich, E., 2006, "Torsional Electromechanical Quantum Oscillations in Carbon Nanotubes," Nature, 1, pp. 36-41. [19] Gurtin, M. E., and Murdoch, A. I., 1975, "A Continuum Theory of Elastic Material Surfaces," Archive for Rational Mechanics and Analysis, 57(4), pp. 291-323. [20] Cuenot, S., Frétigny, C., Demoustier-Champagne, S., and Nysten, B., 2004, "Surface Tension Effect on the Mechanical Properties of Nanomaterials Measured by Atomic Force Microscopy," Physical Review B, 69, pp. 165410.1-165410.5. [21] Dingreville, R., Qu, J., and Cherkaoui, M., 2005, "Surface Free Energy and its Effect on the Elastic Behavior of Nano-Sized Particles, Wires and Films," Journal of the Mechanics and Physics of Solids, 53(8), pp. 1827-1854. [22] McFarland, A. W., Poggi, M. A., Doyle, M. J., Bottomley, L. A., and Colton, J. S., 2005, "Influence of Surface Stress on the Resonance Behavior of Microcantilevers," Applied Physics Letters, 87, pp. 053505. [23] Wang, G. F., Feng, X. Q., and Yu, S. W., 2007, "Surface Buckling of a Bending Microbeam due to Surface Elasticity," Europhysics Letters, 77(4), pp. 44002. [24] He, L., Lim, C., and Wu, B., 2004, "A Continuum Model for Size-Dependent Deformation of Elastic Films of Nano-Scale Thickness," International Journal of Solids and Structures, 41(3-4), pp. 847-857. 94 [25] Baeurle, S. A., Usami, T., and Gusev, A. A., 2006, "A New Multiscale Modeling Approach for the Prediction of Mechanical Properties of Polymer-Based Nanomaterials," Polymer, 47(26), pp. 8604-8617. [26] Xiao, S. P., and Belytschko, T., 2004, "A Bridging Domain Method for Coupling Continua with Molecular Dynamics," Computer Methods in Applied Mechanics and Engineering, 193(17-20), pp. 1645-1669. [27] Qian, D., Wagner, G. J., and Liu, W. K., 2004, "A Multiscale Projection Method for the Analysis of Carbon Nanotubes," Computer Methods in Applied Mechanics and Engineering, 193(17-20), pp. 1603-1632. [28] Xiao, S., and Hou, W., 2004, "Studies of Nanotube-Based Resonant Oscillators through Multiscale Modeling and Simulation," Physical Review B, 75(12), pp. 125414. [29] Eringen, A. C., 1972, "Linear Theory of Nonlocal Elasticity and Dispersion of Plane Waves," International Journal of Engineering Science, 10(5), pp. 425-435. [30] Peddieson, J., Buchanan, G. R., and McNitt, R. P., 2003, "Application of Nonlocal Continuum Models to Nanotechnology," International Journal of Engineering Science, 41(3-5), pp. 305-312. [31] Zhang, Y. Q., Liu, G. R., and Wang, J. S., 2004, "Small-Scale Effects on Buckling of Multiwalled Carbon Nanotubes Under Axial Compression," Physical Review B, 70(20), pp. 205430. [32] Zhang, Y. Q., Liu, G. R., and Han, X., 2006, "Effect of Small Length Scale on Elastic Buckling of Multi-Walled Carbon Nanotubes Under Radial Pressure," Physics Letters A, 349(5), pp. 370-376. [33] Sudak, L. J., 2003, "Column Buckling of Multiwalled Carbon Nanotubes using Nonlocal Continuum Mechanics," Journal of Applied Physics, 94, pp. 7281-7287. [34] Zhang, Y. Q., Liu, G. R., and Xie, X. Y., 2005, "Free Transverse Vibrations of Double-Walled Carbon Nanotubes using a Theory of Nonlocal Elasticity," Physical Review B, 71(19), pp. 195404. [35] Xu, M., 2006, "Free Transverse Vibrations of Nano-to-Micron Scale Beams," Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 462(2074), pp. 2977-2995. [36] Wang, Q., 2005, "Wave Propagation in Carbon Nanotubes Via Nonlocal Continuum Mechanics," Journal of Applied Physics, 98, pp. 124301. 95 [37] Wang, Q., Zhou, G., and Lin, K., 2006, "Scale Effect on Wave Propagation of Double-Walled Carbon Nanotubes," International Journal of Solids and Structures, 43(20), pp. 6071-6084. [38] Hu, Y. G., Liew, K., Wang, Q., He, X., and Yakobson, B., 2008, "Nonlocal Shell Model for Elastic Wave Propagation in Single-and Double-Walled Carbon Nanotubes," Journal of the Mechanics and Physics of Solids, 56(12), pp. 3475-3485. [39] Nanorex Inc., 2005, NanoHive-1 v.1.2.0-b1, www.nanorex.com [40] Eringen, A. C., 1977, "Screw Dislocation in Non-Local Elasticity," Journal of Physics D-Applied Physics, 10, pp. 671-678. [41] Eringen, A. C., 1983, "On Differential Equations of Nonlocal Elasticity and Solutions of Screw Dislocation and Surface Waves," Journal of Applied Physics, 54(9), pp. 4703-4710. [42] Povstenko, Y. Z., 1999, "The Nonlocal Theory of Elasticity and its Applications to the Description of Defects in Solid Bodies," Journal of Mathematical Sciences, 97(1), pp. 3840-3845. [43] Timoshenko, S., and Gere, J.M., 1961, Theory of elastic stability, McGraw-Hill, New York, USA, pp. 452-509. [44] Donnell, L.H., 1976, Beams, plates and shells, McGraw-Hill, New York, USA, pp. 377-445. [45] Alder, B. J., and Wainwright, T. E., 1957, "Phase Transition for a Hard Sphere System," The Journal of Chemical Physics, 27, pp. 1208-1209. [46] Ercolessi, F., 1997, "A Molecular Dynamics Primer," Spring College in Computational Physics, ICTP, Trieste, , pp. 24-25. [47] Allinger, N.L., and Burkert, U., 1992, Accurate Molecular Structures: Their Determination and Importance, Oxford University Press, USA, pp. 336-355. [48] Lennard-Jones, J. E., 1924, "On the Determination of Molecular Fields II. from the Equation of State of a Gas," Proceedings of the Royal Society of London, The Royal Society, 106, pp. 463-477. [49] Tersoff, J., 1989, "Modeling Solid-State Chemistry: Interatomic Potentials for Multicomponent Systems," Physical Review B, 39(8), pp. 5566-5568. [50] Brenner, D. W., 1990, "Empirical Potential for Hydrocarbons for use in Simulating the Chemical Vapor Deposition of Diamond Films," Physical Review B, 42(15), pp. 9458-9471. 96 [51] Dumitrica, T., Hua, M., and Yakobson, B. I., 2006, "Symmetry-, Time-, and Temperature-Dependent Strength of Carbon Nanotubes," National Academy of Sciences, 103, pp. 6105-6109. [52] Berendsen, H. J. C., Postma, J. P. M., Van Gunsteren, W. F., DiNola, A., and Haak, J. R., 1984, "Molecular Dynamics with Coupling to an External Bath," The Journal of Chemical Physics, 81(8), pp. 3684-3690. [53] Mylvaganam, K., and Zhang, L., 2004, "Important Issues in a Molecular Dynamics Simulation for Characterising the Mechanical Properties of Carbon Nanotubes," Carbon, 42(10), pp. 2025-2032. [54] Dumitrică, T., and Yakobson, B. I., 2004, "Strain-Rate and Temperature Dependent Plastic Yield in Carbon Nanotubes from Ab Initio Calculations," Applied Physics Letters, 84(15), pp. 2775-2777. [55] Stuart, S. J., Tutein, A. B., and Harrison, J. A., 2000, "A Reactive Potential for Hydrocarbons with Intermolecular Interactions," The Journal of Chemical Physics, 112, pp. 6472-6486. [56] Liu, J. Z., Zheng, Q. S., Wang, L. F., and Jiang, Q., 2005, "Mechanical Properties of Single-Walled Carbon Nanotube Bundles as Bulk Materials," Journal of the Mechanics and Physics of Solids, 53(1), pp. 123-142. [57] Scarpa, F., and Adhikari, S., 2008, "A Mechanical Equivalence for Poisson's Ratio and Thickness of C-C Bonds in Single Wall Carbon Nanotubes," Journal of Physics- London-D Applied Physics, 41, pp. 085306. [58] Yakobson, B., Brabec, C., and Bernholc, J., 1996, "Nanomechanics of Carbon Tubes: Instabilities Beyond Linear Response," Physical Review Letters, 76(14), pp. 2511-2514. [59] Srivastava, D., Wei, C., and Cho, K., 2003, "Nanomechanics of Carbon Nanotubes and Composites," Applied Mechanics Reviews, 56(2), pp. 215. [60] Kudin, K. N., Scuseria, G. E., and Yakobson, B. I., 2001, "C~ 2F, BN, and C Nanoshell Elasticity from Ab Initio Computations," Physical Review B, 64(23), pp. 235406. [61] Lu, P., Lee, H. P., Lu, C., and Zhang, P. Q., 2007, "Application of Nonlocal Beam Models for Carbon Nanotubes," International Journal of Solids and Structures, 44(16), pp. 5289-5300. [62] Wang, Q., and Varadan, V. K., 2007, "Application of Nonlocal Elastic Shell Theory in Wave Propagation Analysis of Carbon Nanotubes," Smart Materials and Structures, 16(1), pp. 178-190. 97 [63] Wang, C. M., Zhang, Y. Y., and He, X. Q., 2007, "Vibration of Nonlocal Timoshenko Beams," Nanotechnology, 18(10), pp. 105401. [64] Heireche, H., Tounsi, A., and Benzair, A., 2008, "Scale Effect on Wave Propagation of Double-Walled Carbon Nanotubes with Initial Axial Loading," Nanotechnology, 19(18), pp. 185703. [65] Kumar, D., Heinrich, C., and Waas, A. M., 2008, "Buckling Analysis of Carbon Nanotubes Modeled using Nonlocal Continuum Theories," Journal of Applied Physics, 103, pp. 073521. [66] Lu, P., Lee, H. P., Lu, C., and Zhang, P. Q., 2006, "Dynamic Properties of Flexural Beams using a Nonlocal Elasticity Model," Journal of Applied Physics, 99, pp. 073510. [67] Sears, A., and Batra, R. C., 2004, "Macroscopic Properties of Carbon Nanotubes from Molecular-Mechanics Simulations," Physical Review B, 69(23), pp. 235406. [68] Huang, Y., Wu, J., and Hwang, K. C., 2006, "Thickness of Graphene and Single- Wall Carbon Nanotubes," Physical Review B, 74(24), pp. 245413. 98 Appendix This appendix presents a brief overview of the AIREBO potential. For full details, the reader is encouraged to consult the original paper by Stuart et al. [55]. The AIREBO potential is represented by the sum of all the pairwise interactions in the system, which includes the covalent bonding interactions (REBO potential) as well as Lennard-Jones terms and torsional interactions: , , , 1 2 AIREBO REBO LJ tors ij ij kijl i j i k i j l i j k E E E E ≠ ≠ ≠ ⎡ ⎤= + +⎢ ⎥⎢ ⎥⎣ ⎦ ∑∑ ∑ ∑ . (A.1) The REBO interactions are expressed by: ( ) ( )REBO R Aij ij ij ij ij ijE V r b V r= + , (A.2) where RijV and A ijV are repulsive and attractive pair-wise potentials, rij is the inter-atomic distance and bij is the many-body (or bond-order) term. The repulsive term is: ( ) ( ) 1 ij ijrijRij ij ij ij ij ij Q V r w r A e r α−⎡ ⎤= +⎢ ⎥⎢ ⎥⎣ ⎦ , (A.3) where Qij, Aij and αij depend on atom types i and j. The values of these parameters can be found in reference [55]. wij is a bond-weighting factor: ( ) ( ( ))ij ij c ijw r S t r′= . (A.4) The bond-weighting factor switches off the REBO interactions when the atom pairs exceed typical bonding distances. The switching function is expressed as: 99 [ ]1( ) ( ) ( ) (1 ) 1 cos( ) 2 S t t t t tπ′ = Θ − + Θ Θ − + , (A.5) and the switching region for each type of bond is given by a scaling function: ( ) min ij ij c ij max min ij ij r r t r r r −= − . (A.6) The attractive term in equation (A.2) is given by: ( )3 ( ) 1 ( ) ( ) n ijij rA n ij ij ij ij ij n V r w r B e β− = = − ∑ . (A.7) The bij term in equation (A.2) specifies the bond-order for the interaction between atoms i and j: 1 2 rc dh ij ij ji ij ijb p p σπ σπ π π⎡ ⎤= + + +⎣ ⎦ . (A.8) The bij term is larger for stronger bonds and each of the terms that contribute to bij is a many-body term that depends on the bonding environment surrounding atoms i and j. The principal contributor is the covalent bonding interaction: 1/2 , 1 ( ) (cos ) jikij ik ik i ijk ij k i j p w r g e Pλσπ θ − = ⎡ ⎤= + +⎢ ⎥⎢ ⎥⎣ ⎦ ∑ , (A.9) where θjik is the bond angle between rji and rki vectors. gi is a penalty function that imposes a cost on bonds that are too close to one another and depends on the local coordination number, defined as the sum of carbon-only and hydrogen-only coordination numbers: 100 C H ij ij ijN N N= + , (A.10) where: ( ) ( )Cij kC ik ik jC ij ij k i N w r w rδ δ ≠ ⎛ ⎞= −⎜ ⎟⎝ ⎠∑ (A.11) is the carbon-only coordination number and δij is the Kronocker delta. The penalty function is thus defined as: (1) (2) (1)(cos ) (cos ) ( ( )) (cos ) (cos )C jik jik N ij jik jikC C Cg g S t N g gθ θ θ θ⎡ ⎤′= + −⎣ ⎦ . (A.12) tN is a scaling function given by: ( ) min ij ij N ij max min ij ij N N t N N N −= − . (A.13) The remaining terms in equation (A.9) are small correction factors which are described in reference [55]. Looking back at equation (A.8), the term rcijπ includes the contributions to the bond order from the radical and conjugation effects. The local measure of conjugation in the bond i-j is: 2 2 , , 1 ( ) ( ( )) ( ) ( ( ))conjij kC ik ik conj ki lC jl jl conj lj k i j l i j N w r S t N w r S t Nδ δ ≠ ≠ ⎡ ⎤ ⎡ ⎤′ ′= + +⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎣ ⎦ ⎣ ⎦ ∑ ∑ , (A.14) where: 101 ( ) min conj max min N Nt N N N −= − . (A.15) The final term contributing to the bond order bij is dhijπ , which imposes a penalty for rotation around multiple bonds, a detailed description of which is given in reference [55]. The Lennard-Jones contribution to the pair energy is expressed by: *( ( )) ( ( )) ( ) 1 ( ( )) ( )LJ LJ LJij r ij b ij ij ij ij r ij ij ij ijE S t r S t b C V r S t r C V r⎡ ⎤= + −⎣ ⎦ , (A.16) where ( )LJij ijV r is the traditional Lennard-Jones term: 12 6 ( ) 4 ij ijLJij ij ij ij ij V r r r σ σε ⎡ ⎤⎛ ⎞ ⎛ ⎞⎢ ⎥= −⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦ . (A.17) The switching function S(t) is defined as: 21( ) ( ) ( ) (1 ) 1 (3 2 ) 2 LJ ijE S t t t t t t⎡ ⎤= = Θ − + Θ Θ − − −⎣ ⎦ . (A.18) tr is the scaling function defined as: ( ) LJ min ij ij r ij LJ max LJ min ij ij r r t r r r −= − . (A.19) At intermolecular distances, the LJ interaction is included only if there is no significant bonding interaction between the two atoms and if the atoms i and j are not connected by two or fewer intermediate atoms, as specified by the tb and Cij switches, respectively: 102 ( ) min ij ij b ij max min ij ij b b t b b b −= − , (A.20) ( ), ( ) ( ), 1 max ( ) ( ) ( ), , ij ij ik ik kj kj ij ik ik ik ik ik ik w r w r w r k C w r w r w r k l ∀⎧ ⎫⎪ ⎪= − ⎨ ⎬∀⎪ ⎪⎩ ⎭ . (A.21) Finally, the torsional potential of the AIREBO potential is due to the dihedral angle determined by atoms i,j,k and l: ( ) ( ) ( ) ( )tors torskijl ki ki ij ij jl jl kijlE w r w r w r V w= , (A.22) where: 10256 1cos ( / 2) 405 10 tors kijl kijl kijl kijlV wε ε= − . (A.23)