Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 50, 4, pp. 975-986, Warsaw 2012 50th Anniversary of JTAM BUCKLING ANALYSIS OF SHORT CARBON NANOTUBES BASED ON A NOVEL TIMOSHENKO BEAM MODEL Reza Hosseini-Ara, Hamid R. Mirdamadi, Hasan Khademyzadeh Isfahan University of Technology, Department of Mechanical Engineering, Isfahan, Iran e-mail: r.hosseiniara@me.iut.ac.ir In this paper, we present a novelmethod to investigate the buckling behavior of short clam- ped carbon nanotubes (CNTs) with small-scale effects. Based on the nonlocal Timoshenko beam kinematics, the strain gradient theory and variational methods, the higher-order go- verning equation and its correspondingboundary conditions are derived,which are often not considered. Then, we solve the governing differential equation and determine exact critical buckling loads using a linear polynomial plus trigonometric functions different from the pu- rely trigonometric series. We also investigate the influences of the scale coefficients, aspect ratio and transverse shear deformation on the buckling of short clamped CNTs. Moreover, we compare the critical strains with the results obtained from the Sanders shell theory and validate themwithmolecular dynamic simulationswhich are found to be in good agreement. The results show that unlike the other beam theories, this model can capture correctly the small-scale effects on buckling strains of short CNTs for the shell-type buckling. Key words: short carbon nanotubes, buckling, nonlocal elasticity, Timoshenko beam theory 1. Introduction Since discovering carbon nanotubes (CNTs) by Iijima (1991) more than a decade ago, scientists have been exploring possible uses for carbonnanotubes,which exhibit superior electrical, chemi- cal andmechanical properties. Owing to these properties, CNTs can be applied as nano-probes, nano-needles, reinforcing fibers in composites, nano-actuators, nano-vessels for hydrogen storage and gene delivery systems and nanoscale electronic devices in nano-electromechanical systems (Zhang et al., 2009; Wang et al., 2010). Meanwhile, the application of short CNTs has prompted significant effort to reduce the size of these nanoscale devices. Seidel et al., (2005) showed that short CNTs, with lengths less than 20 nanometers, are useful inmolecular electronics andCNTfield-effect transistors (CNTFETs). In addition, they are smaller than many large proteins in the bloodstream, so tubes of that length could find uses as biomedical sensors. As a result, the importance of carbon nanotubes is realized and both theoretical and experi- mentalworks are carried out (Muc, 2011;Yakobson et al., 1996).Oneproperty in particular that has been extensively studied is the buckling of single-walled carbon nanotubes (SWCNTs) un- der axial compression. In fact, nanotubes are highly susceptible to buckling under compression, which is a structural instability. Once the buckling of CNTs occurs, the load-carrying capabi- lity would suddenly reduce and lead to possible catastrophic failure of the nanotubes, which can significantly limit the loading strengths of the probing tips and compressive strengths of nanocomposites. Even the physical properties such as conductance of a carbon nanotube can be influenced by the occurrence of buckling (Postma et al. 2001). Hence, it is crucial to understand themechanism of nanotube buckling and even predict the onset of buckling in order to improve the nanotube applications. 976 R. Hosseini-Ara et al. To this end, continuummechanical modelingmethods andMolecular dynamics (MD) simu- lations therefore play important roles in investigating the key characters of the structural and mechanical properties of CNTs.Yakobson et al. (1996) conducted some of the early studies that explored the buckling behavior of SWCNTs under various conditions including bending, com- pression, and twisting. They proposed that classical mechanics can be employed to predict the buckling strain of SWCNTsunder strainbaseduponcertain approximations.However, at atomic length scale, thematerial microstructure becomes increasingly important and small-scale effects cannot be ignored. In order to improve the CNT constitutive relations, many authors adopted Eringen’s equations of nonlocal elasticity introduced byEringen andEdelen (1972) andEringen (1972, 1976, 1983, 2002), and incorporated them into several continuummodels (Feliciano et al., 2011; Kumar et al., 2008; Lu et al., 2006; Peddieson et al., 2003; Reddy, 2007; Reddy and Pang, 2008; Sudak, 2003; Wang B. et al., 2010; Wang C.M., et al., 2007; Wang Q. et al., 2006; Wang Q., 2005; Wang Q. and Wang C.M., 2007; Wang Y.Z. et al., 2010; Yang et al., 2010; Zhang et al., 2004). The use of nonlocal elasticity relations in either Euler-Bernoulli or Timoshenko beam models were shown to be accurate for the buckling strains of the long CNTs (Silvestre et al., 2011; Zhang et al., 2009). In addition, for short SWCNTs with large diameters, the nonlocal shell model with the appropriate small length scale parameter can provide critical strains that are in good agreement withMD results (Silvestre et al., 2011). However, for short SWCNTswith small diameters,more work has to be done to refine the nonlocal beamand shellmodels for better prediction of critical strains (Silvestre et al., 2011). In this paper, we investigate the use of the more refined Timoshenko beam theory for mo- deling the buckling behavior of CNTswith small aspect ratios. The nonlocal Timoshenko beam model results are compared with the Sanders thin shell model results and validated against MD simulations results. It will be shown herein that unlike the other beammodels, the current Timoshenko beam model can correctly reproduce the buckling strains of short CNTs that are length dependent, and these results are relatively close to those predicted byMD simulations. 2. Nonlocal Timoshenko beam theory The nonlocal elasticity model was first presented by Eringen (1972). According to this model, the stress at a reference point in the body is dependent not only on the strain state at that point, but also on the strain state at all of the points throughout the body. The constitutive equation of the nonlocal elasticity can be written as follows [1− (e0a)2∇2]σij =Cijklεkl (2.1) where Cijkl is the elastic modulus tensor of the classical isotropic elasticity; and σij and εkl are the stress and strain tensors, respectively. In addition, e0 is a nondimensionalmaterial constant, determined by experiments, and a is an internal characteristic length (e.g., a lattice parameter, granular distance). Therefore, e0a is a constant parameter that showing the small-scale effect in nano-structures. The assumed displacement field of the Timoshenko beam kinematics is u1(x,y,z,t) =u(x,t)+zφ(x,t) u2(x,y,z,t) = 0 u3(x,y,z,t) =w(x,t) (2.2) where φdenotes the rotation of the cross section at the point x about the y-axis. The remaining nonzero axial and transverse shear strains are given by εxx = ∂u ∂x +z ∂φ ∂x γ≡ 2εxz = ∂w ∂x +φ (2.3) Buckling analysis of short carbon nanotubes... 977 Using Eq. (2.1), the nonlocal stress tensor components are σxx− (e0a)2 ∂2σxx ∂x2 =Eεxx σxz− (e0a)2 ∂2σxz ∂x2 =KSGγ (2.4) The nonlocal stress resultants of the axial force, shear, and bending moment are derived from the above equations, respectively NNL− (e0a)2 ∂2NNL ∂x2 =EA ∂u ∂x QNL− (e0a)2 ∂2QNL ∂x2 =KSGAγ MNL− (e0a)2 ∂2MNL ∂x2 =EI ∂φ ∂x (2.5) where E, G, A and I are Young’s modulus, shear modulus, cross-sectional area of beam and area moment of inertia of the beam cross section, respectively. Additionally, KS denotes the shear correction factor, defined by Q=KS ∫ A σxz dA (2.6) This factor corrects the assumption of constant shear strain on the cross section of the beam in the Timoshenkomodel, depending on thematerial and geometry of the cross section. 3. Strain gradient elasticity Solving Eq. (2.4), the nonlocal axial and shear stresses as a function of the displacement field can be determined as follows σ(x,z) =E [ ε(x,z)+(e0a) 2∂ 2ε(x,z) ∂x2 +(e0a) 4∂ 4ε(x,z) ∂x4 + . . . ] τ(x)=KSG [ γ(x)+(e0a) 2d 2γ(x) dx2 +(e0a) 4d 4γ(x) dx4 + . . . ] (3.1) Assuming (e0a) 2 ≪ 1, and neglecting the higher powers of the nonlocal parameter, (e0a)4, (e0a) 6, etc., the solution could be simplified to σ(x,z) =E [ ε(x,z)+(e0a) 2∂ 2ε(x,z) ∂x2 ] τ(x)=KSG [ γ(x)+(e0a) 2d 2γ(x) dx2 ] (3.2) In fact, Eqs. (3.2) can be thought of as constituting a strain gradient form of the nonlocal beam model (Peddieson et al., 2003).Considering the straingradient approach, for aTimoshenkobeam subjected to an external compressive and conservative force field N0, and laterally distributed load p(x), the total potential energy Π, is given by Kumar et al. (2008) is generalized in the presence of shear deformation effects as follows Π = ∫ V [Eε2(x,z) 2 − (e0a)2 E 2 (∂ε(x,z) ∂x )2] dv+ ∫ V [KSGγ 2(x) 2 − (e0a)2 KSG 2 (dγ(x) dx )2] dv + L ∫ 0 ( N0 du dx ) dx− L ∫ 0 pw(x) dx− 1 2 L ∫ 0 N0 (dw dx )2 dx (3.3) where N0 is an external and axial compressive load. Thus, the second term is added to the original equation for capturing the shear deformation effects in Timoshenko beam theory. Fur- thermore, Chang et al. (2002) prove the original form of this equation for strain gradient theory 978 R. Hosseini-Ara et al. without higher-order stress. They used the characteristic size coefficient (d2/6) instead of the nonlocal parameter and derived the potential energy density using integration by parts. The last three terms of Eq. (3.3) are also the work done by the axial load, lateral load and VonKarman effect, respectively. Substituting Eq. (2.3) into Eq. (3.3) and integrating over the cross-sectional area, the follo- wing expression is obtained for Π Π = 1 2 L ∫ 0 [ EA (du dx )2 +EI (dφ dx )2 +KSGA (dw dx +φ )2] dx − 1 2 (e0a) 2 L ∫ 0 [ EA (d2u dx2 )2 +EI (d2φ dx2 )2 +KSGA (d2w dx2 + dφ dx )2] dx + L ∫ 0 ( N0 du dx ) dx− L ∫ 0 pw(x) dx− 1 2 L ∫ 0 N0 (dw dx )2 dx (3.4) 4. Governing differential equations and boundary conditions The classical axial force NCL acting on the beam cross-section is defined as NCL =EA du dx (4.1) Using the above expression for NCL and ignoring the laterally distributed loads p for buckling analysis, the variation of Eq. (3.4) with respect to u(x) and equating to zero, it can be written as δu(Π)= L ∫ 0 [ NCL dδu dx − (e0a)2 (dNCL dx )d2δu dx2 +N0 dδu dx ] dx=0 (4.2) Integrating by parts, we get δu(Π)= [NCLδu] x=L x=0 − L ∫ 0 [dNCL dx δu ] dx− [ (e0a) 2 (dNCL dx )dδu dx ]x=L x=0 + L ∫ 0 [ (e0a) 2 (d2NCL dx2 )dδu dx ] dx+ [ N0δu ]x=L x=0 − L ∫ 0 [dN0 dx δu ] dx=0 (4.3) Integrating by parts again, we obtain δu(Π)= [NCLδu] x=L x=0 − L ∫ 0 [d,NCL dx δu ] dx− [ (e0a) 2 (dNCL dx )dδu dx ]x=L x=0 + [ (e0a) 2 (d2NCL dx2 ) δu ]x=L x=0 − L ∫ 0 [ (e0a) 2 (d3NCL dx3 ) δu ] dx + [ N0δu ]x=L x=0 − L ∫ 0 [dN0 dx δu ] dx=0 (4.4) Buckling analysis of short carbon nanotubes... 979 Assuming the strain gradient approach in Eq. (3.2)1 for the nonlocal axial force NNL and by some simplification, we obtain δu(Π)= [ NCL+(e0a) 2 (d2NCL dx2 ) +N0 ] δu ∣ ∣ ∣ ∣ x=L x=0 − L ∫ 0 [dNCL dx +(e0a) 2 (d3NCL dx3 ) + dN0 dx ] δu dx− [ (e0a) 2 (dNCL dx )dδu dx ]x=L x=0 =(NNL+N0)δu ∣ ∣ ∣ x=L x=0 − [ (e0a) 2 (dNCL dx )dδu dx ]x=L x=0 − L ∫ 0 [dNNL dx + dN0 dx ] δu dx=0 (4.5) Thus, the governing equation and boundary conditions in the x direction are derived fromEqs. (4.2)-(4.5), as follows d dx (NNL+N0)=0 (NNL+N0)δu ∣ ∣ ∣ x=L x=0 =0 − (e0a)2 (dNCL dx ) δ du dx ∣ ∣ ∣ ∣ x=L x=0 =0 (4.6) Performing variation with respect to w(x) for Eq. (3.4) and equating to zero, it gives δw(Π)= L ∫ 0 [ KSGA (dw dx +φ )dδw dx − (e0a)2KSGA (d2w dx2 + dφ dx )d2δw dx2 ] dx − L ∫ 0 ( N0 dw dx dδw dx ) dx=0 (4.7) Integrating by parts, we obtain the governing equation and boundary conditions for w as dQNL dx =N0 d2w dx2 [ QNL−N0 dw dx ] δw ∣ ∣ ∣ ∣ L 0 =0 [ −(e0a)2KSGA (d2w dx2 + dφ dx )]dδw dx ∣ ∣ ∣ ∣ L 0 =0 (4.8) In the same way, applying the variational operator to φ(x) for Eq. (3.4) and equating to zero, we obtain δφ(Π)= L ∫ 0 [( EI dφ dx )dδφ dx − ( (e0a) 2EI d2φ dx2 )d2δφ dx2 +KSGA (dw dx +φ ) δφ ] dx − L ∫ 0 [ (e0a) 2KSGA (d2w dx2 + dφ dx )dδφ dx ] dx=0 (4.9) Using integration by parts, the governing equation is given by dMNL dx =QNL (4.10) and the following boundary conditions are derived [ MNL− (e0a)2KSGA (d2w dx2 + dφ dx )] δφ ∣ ∣ ∣ ∣ L 0 =0 [ −(e0a)2EI d2φ dx2 ]dδφ dx ∣ ∣ ∣ ∣ L 0 =0 (4.11) 980 R. Hosseini-Ara et al. Substituting the nonlocal shear force and bending moment defined in Eqs. (2.5)2,3 into the governingEqs. (4.8)1 and (4.10), andomitting the similar terms fromboth sides of the equations, we obtain KSGA (d2w dx2 + dφ dx ) +(e0a) 2N0 d4w dx4 =N0 d2w dx2 EI d2φ dx2 =KSGA (dw dx +φ ) (4.12) Solving Eq. (4.12)1 for φ gives dφ dx = 1 KSGA [ N0 d2w dx2 − (e0a)2N0 d4w dx4 ] − d2w dx2 (4.13) By differentiating Eq. (4.12)2 and inserting Eq. (4.13) in Eq. (4.12)2, the transverse equilibrium equation in terms of the lateral displacement for an axially loaded beam using a nonlocal strain gradient theory is obtained as [EI(e0a) 2N0 KSGA ]d6w dx6 + [ EI− EI KSGA N0− (e0a)2N0 ]d4w dx4 +N0 d2w dx2 =0 (4.14) where N0 is an external axial compressive load. This equation is similar to that obtained by Reddy and Pang (2008) for buckling of the nonlocal Timoshenko beam using Hamilton’s Prin- ciple. In addition, for solving the above equation, six boundary conditions are required (three for each end) but eight boundary conditions appear in Eqs. (4.8)2,3 and (4.11). It means that there is one additional boundary condition for each end. So, the main objective is to select three independent boundary conditions which can satisfy all four boundary conditions for each end. In the next part, the boundary conditions for various beam supports are obtained. Thenondimensional formofEq. (4.14) using length of the beam L, as a nondimensionalizing parameter can be rewritten as (L4µΩ) d6w dx6 +L2 ( 1 π2r −Ω−µ )d4w dx4 + d2w dx2 =0 (4.15) where Ω and µ are the nondimensional forms of shear deformation and nonlocal parameters, respectively, and r is the ratio of the critical buckling loads as follows Ω= EI KSGAL2 µ= (e0a L )2 r= NNLcr NLcr (4.16) where NNLcr is obtained from solving the Eq. (4.14) and N L cr is that given by classic Euler columns for simply supported end conditions. We may simply switch to nonlocal Euler-Bernoulli beam model by ignoring the shear de- formation terms. Additionally, the local Timoshenko beam model is obtained by letting the nonlocal parameter to be zero and by setting the shear deformation and nonlocal parameters to zero, the local Euler-Bernoulli beammodel appears. 5. Buckling solutions Here we consider analytical solutions based on a linear polynomial plus trigonometric functions different from the purely trigonometric series, for nonlocal Timoshenko beams under a constant axial compressive load, using the buckling equation obtained in Eq. (4.14) for the clamped end Buckling analysis of short carbon nanotubes... 981 conditions. This sixth order equation exhibits different solutions depending on the ratio r, and the nonlocal and shear deformation parameters µ and Ω. The discriminant of the characteristic equation corresponding to differential Eq. (4.14) is defined as follows ∆= [EI N0 − EI KSGA − (e0a)2 ]2 −4 [EI(e0a) 2 KSGA ] =L4 [( 1 π2r −Ω−µ )2 −4µΩ ] (5.1) If ∆> 0, one of the solutions is w(x) = c1+ c2x+c3 sin(Px)+ c4cos(Px)+ c5 sin(Qx)+ c6cos(Qx) (5.2) where c1,c2, . . . ,c6 are constants of integration and determined by six boundary conditions. P and Q are given by P = √ √ √ √ √ EI N0 − EI KSGA − (e0a)2− √ ∆ 2 [ EI(e0a)2 KSGA ] Q= √ √ √ √ √ EI N0 − EI KSGA − (e0a)2+ √ ∆ 2 [ EI(e0a)2 KSGA ] (5.3) If ∆< 0, the other solution is defined as w(x) = c1+ c2x+c3e Rx cos(Sx)+ c4e Rx sin(Sx)+ c5e −Rx cos(Sx)+ c6e −Rx sin(Sx) (5.4) where R and S are R= 1 2      KSGA EI(e0a)2  2+ EI KSGA +(e0a) 2− EI N0 e0a √ EI KSGA   2      1 4 S = 1 2      KSGA EI(e0a) 2  2− EI KSGA +(e0a) 2− EI N0 e0a √ EI KSGA   2      1 4 (5.5) The first solution in Eq. (5.2) is found to be valid depending on the sign of ∆. Therefore, the second solution is not used in this research and is only stated for completeness. 5.1. Clamped beams Regarding theclassical continuummechanics for clampedboundaryconditions, thedeflection and rotation of the cross-section are zero at each boundary. Thus, we use these classical as well as the newly derived boundary conditions in Eqs. (4.8)2,3 and (4.11) to derive the following boundary conditions. To this end, we satisfy all four boundary conditions, as defined in Eqs. (4.8)2,3 and (4.11) for each end w=0 dw dx =0 [N0(e0a) 2] d5w dx5 +(KSGA−N0) d3w dx3 =0 (5.6) Substituting theseboundaryconditions inEq. (5.2) and setting thedeterminant of the coefficient matrix to be zero, the critical buckling load is derived as Ncritical = 4π2L2EIKSGA L4KSGA+4π2L2EI+4π2L2(e0a)2KSGA+16π4EI(e0a)2 (5.7) The critical buckling load using the dimensionless parameters µ and Ω is Ncritical = 4π2EI L2 [ 1 1+4π2(µ+Ω)+16π4(µΩ) ] (5.8) This equation may be transformed into the nonlocal Euler-Bernoulli theory, for the shear para- meter set to zero, Ω =0, or classical Timoshenko beam theory, for the nonlocal parameter set to zero, µ = 0, or classical Euler-Bernoulli beam theory, for both the shear and the nonlocal parameters set to zero, Ω=µ=0. 982 R. Hosseini-Ara et al. 6. Numerical results 6.1. Comparison of critical buckling loads for beam theories In this section,weconsidernumerical solutions forCNTsmodeledasnanobeamswith circular cross sections. The numerical results are presented in the form of graphs and tables for different types of end conditions, using the following effective properties of carbonnanotubes (Reddy and Pang, 2008) E =1000GPa G=420GPa d=1nm KS =0.877 I = πd4 64 =0.049(nm)4 A= πd2 4 =0.785(nm)2 (6.1) A plot of the critical buckling loads for the nonlocal clamped Timoshenko beam with different values of shear deformation and nonlocal parameters is presented in Fig. 1. Fig. 1. Plot of the ratios of buckling loads for different values of µ and Ω As illustrated inFig. 1, the solid lines for Ω=0denoteEuler-Bernoulli beamswhich are the upper bound solutions. By increasing the shear deformation parameter, Ω, the critical buckling loads decrease. The effect of shear deformation is quantified for different boundary conditions. This effect is negligible for L/d ratios more than 20 (or Ω less than 0.0005), but significant by increasing the Ω, for L/d ratios less than 20. As we are discussing about the short CNTs (i.e., 2