Stability Analysis of Nonlocal Elastic Columns with Initial Imperfection

Investigated herein is the postbuckling behavior of an initially imperfect nonlocal elastic column, which is simply supported at one end and subjected to an axial force at the other movable end. The governing nonlinear differential equation of the axially loaded nonlocal elastic column experiencing large deflection is first establishedwithin the framework of Eringen’s nonlocal elasticity theory in order to embrace the size effect. Its semianalytical solutions by the virtue of homotopy perturbation method, as well as the successive approximation algorithm, are determined in an explicit form, through which the postbuckling equilibrium loads in terms of the end rotation angle and the deformed configuration of the column at this end rotation are predicted. By comparing the degenerated results with the exact solutions available in the literature, the validity and accuracy of the proposed methods are numerically substantiated. The size effect, as well as the initial imperfection, on the buckled configuration and the postbuckling equilibrium path is also thoroughly discussed through parametric studies.


Introduction
The discovery of carbon nanotubes ushers a new era in the nano world [1].The enhanced features of nanomaterial in mechanical, electrical, optical, and chemical properties over traditional ones open up many application opportunities in cutting-edge fields that range from sensing, and communications to energy harvesting.However, when the size of the structure scales down to nanodomains, an issue of considerable importance, namely, the size effect, arises and becomes prominent [2].The size effect may greatly alter the macroscopic properties, and moreover, it calls the applicability of classical continuum models into question.Various modified continuum theories (such as couple stress theory, strain gradient elasticity theory, and modified couple stress theory) have thus been proposed to account for the size effect in the micro/nanoscale structures.Among them, the nonlocal elasticity theory, pioneered by Eringen [3], has been widely accepted and attracted an ever growing attention in recent years.This is because Eringen's nonlocal continuum-based models are mathematically easy to tackle and physically reasonable from the atomistic viewpoint of lattice dynamics and molecular dynamics simulations [4].This theory abandons the classical assumption of locality and admits that the stress state depends not only on the strain at that point but on the strains of every point in the body.In this way, information concerning about the long-range forces between atoms is incorporated into the theory, and consequently, the internal size scale is represented in the constitutive equations simply as a material parameter.
So far, there have been considerable studies on size effects in problems of bending, vibration, buckling, and wave propagation of nanostructures (see [5] and references therein) based on the nonlocal elasticity theory.Among them, buckling is one of the important topics because nanoelements are more susceptible to buckling instability when subjected to compressive loads due to their high aspect ratio.Moreover, the onset of buckling will significantly compromise the structural integrity of nanostructures [6].However, it should be mentioned that most studies assumed small deformations for the nanoelements so that the curvature may be approximated by the second derivative of the transverse deflection.Nevertheless, these nanoelements usually possess excellent flexibility, just as experimental observations on carbon nanotubes have indicated, and they can undertake significant deformations without experiencing apparent plasticity [7][8][9].Thus, the general linearized theories fail to depict this aforementioned behavioral feature in a satisfactory manner [10].In order to predict the postbuckling behavior of nanoelements, both large deflection and size effect must be considered.Hence, a geometrically exact theory is indispensable for postbuckling behavior of nanoelements.
In the light of the mathematical consistency in the way the mechanical model is formulated, hereinafter, we will restrict our attention to the elastica theory [11,12], which entails using the exact expression of curvature.The first study of the deformation of an elastica under loading was probably due to James Bernoulli although it was Euler who established the basic theory [11,13].By taking advantage of Euler's celebrated work, considerable efforts have been devoted to the theory with applications extended from classical fields to statistical mechanics (see [14] and the references therein), biomechanics [15], and ocean engineering [16,17].Within the scope of most of these studies, linear constitutive relations are utilized.Encouraged by its effectiveness in dealing with these geometrical nonlinear problems, the focus is also transferred to other kinds of materials and structures.For example, Haslach Jr. [18] investigated the postbuckling behavior of a column made of a material having a cubic constitutive equation, while Al-Sadder and Shatarat [19] examined the large deflection of a prismatic composite cantilever beam comprising two different nonlinear elastic materials.Within the framework of the elastica theory, the nonlinear buckling response of asymmetric laminated structures composed of an arbitrary number of distinct material layers was treated by Vinogradov and Derrick [20].Kang and Li [21] considered the mechanical behaviors of a cantilevered functionally graded material (FGM) beam that obeys the Ludwick-type law.An interesting "partially composite elastica" problem was studied by Challamel [22], where each subcolumn is modeled with the Euler-Bernoulli beam theory and connected to each other via a linear constitutive law for the interlayer slip.In addition, Frostig [23] revealed the elastica behavior of a sandwich panel with a soft/compliant core that analysis includes both the shear deformation and moderate strain.The large deflection of the fluid-saturated poroelastic beams which are permeable in the axial direction and impermeable in the transverse directions was examined by Li et al. [24].Likewise, an elastica-type dynamic stability of viscoelastic columns with its behavior given in terms of the Boltzmann superposition principle was discussed by Suire and Cederbaum [25].Moreover, the boundary conditions were also extended from the commonly used conditions, namely, simply supported and clamped-free conditions, to clamped-hinged condition [26], frictional end [27], and slippery support [28] conditions.
As for the governing equations of these elastica problems, the elliptical integral functions are often employed to express their exact solutions.For example, by utilizing the elliptical integral transformation, Kimball and Tsai [29] presented the nonlinear deflection solution for a cantilever beam subjected to an arbitrary combination of end forces and moments.Moreover the solution for double curvature bending of variable-arc-length elasticas was proposed by Chucheepsakul et al. [30] in terms of elliptical integrals.However, from an engineering point of view, these exact solutions are not quite intuitive of any physical insight, and one often needs to resort to numerical methods.Besides, the exact treatment remains intractable for most cases because of the intrinsic nonlinearity that resides in the elastica theory.Therefore, another branch of treatment from the numerical point of view is developed simultaneously in conjunction with the availability of high-speed computers.One popular approach is the finite-element method, for example, by which, along with Newton-Rhapson iteration technique, the static analysis of risers under the influence of internal flow and hydrostatic pressure was performed by Chen et al. [16].Another frequently used approach is the shooting method.It converts the nonlinear boundary-value problem into an initial-value problem, which is then determined with or without iterations [31,32].Other newly developed numerical methods include the Lie-group shooting method [33], method of genetic algorithm [34,35], and differential quadrature method [36].In view that an analytical solution is always convenient for parametric studies and captures clearly the physics of the problem, it appears more appealing than the numerical one.Thus, extensive research efforts were also made in this third line of development, which made use of all kinds of asymptotic methods.Among them, the perturbation technique [37][38][39][40][41] is popular in the broad area of buckling problems of structures.For instance, a higher-order perturbation technique was employed by Lacarbonara [39] to determine the postbuckling solutions for nonprismatic nonlinearly elastic rods.In surveying the postbuckling regime for the five classical Euler buckling cases of an extensible elastica, the method of multiple (spatial) scales perturbation was proposed by Mazzilli [40].The analysis of postbuckling, as well as nonlinear bending and nonlinear vibration, for a simply supported Euler-Bernoulli beam resting on a twoparameter elastic foundation was accomplished by Shen [41] with the help of a novel two-step perturbation technique.Another recently developed method, the homotopy analytical method, was used by Wang et al. [42] to solve a cantilever beam problem in order to obtain an explicit solution of the displacement and rotation at the free end.Later, this problem was revisited by Tolou and Herder by means of the Adomian decomposition method [43].The remarkable accuracy of this method was also demonstrated in dealing with the tip loaded cantilever elastica and plastica problems [44].Besides, with the help of a series of admissible functional transformations, Andriotaki et al. [45] furnished an implicit analytical solution for the elastica problem of a cantilever bar due to its own weight.
Although there are extensive studies related to the elastica problems, just as clearly indicated earlier, the scrutinies of large deflections of flexible nanocolumns with size effect included are less frequent.Until recently, efforts toward an elastica type buckling analysis of nanocolumn were conducted.By employing the shooting method, Wang et al. [46] examined the postbuckling problem of cantilevered nanorods/tubes in the absence of shear deformation.Rather than relying on the numerical simulation, Xu [47] independently investigated this problem theoretically with the help of the homotopy perturbation method.For stubby nanocolumns, Xu et al. [48] exploited the Timoshenko hypothesis as a basis for the nonlinear postbuckling description, and in such a way the role of shear deformation in the postbuckling behavior was revealed.By using Pontryagin's maximum principle, the optimal shape of a nonlocal elastic rod clamped at both ends was determined by Atanackovic et al. [49].In these cited articles, the exact curvature expression was used, and the nonlocal elasticity theory was chosen to account for the size effect.
Though meaningful results with allowance for the size effect have been observed, it appears that the postbuckling of nanocolumns still merits some further investigation.In fact, the nanocolumns often suffer from geometric imperfections due to various manufacturing and environmental factors, and the appearance of these initial imperfections in axially loaded structures significantly affects their response in the prebuckling and postbuckling stages [50].On the other hand, initial imperfections may be introduced purposely to take advantage of its beneficial features.For example, the postbuckling behavior of nanocolumn with an initial imperfection to create a very shallow arch can be used as a particular strategy to electronic devices [51].The imperfection amplitude of frames [52] or composite beams (see [53] and references therein) can also be properly manipulated for a desired response, such as maximizing the buckling load.Motivated by those observations, the buckling behavior, especially in the postbuckling stage, of an initially imperfect nonlocal elastic column will be explored in this study.The primary objective will focus on obtaining explicit semianalytical solutions via simple but vigorous approximate techniques.Moreover, it is of interest to evaluate the effect of initial imperfection, as well as the size effect, on the postbuckling behavior of nonlocal elastic columns.
The remainder of this paper is organized as follows.The problem formulation and exact governing equation are furnished in Section 2. Section 3 presents the novel employment of two semianalytical methods, namely, the homotopy perturbation method (HPM) and successive approximate algorithm (SAA), for tackling the challenging geometrically nonlinear problems and obtaining approximate solutions.In Section 4, the main numerical results by the proposed methods are presented and discussed.The concluding remarks are given in Section 5.

Theoretical Formulation
Consider the stability problem of an elastic nonlocal elastic column with a slight geometrical curvature as an imperfection.This inextensible, simply supported column of uniform cross-section  and length  is subjected to a conservative force  at its right movable end as shown in Figure 1(a).The Cartesian (, ) coordinate system is chosen in such a manner that the abscissa axis coincides with the line connecting the two hinged ends, and the coordinate origin is located at its left end.Let  be the arc length reckoned along the column,  the deflection in  direction, and   () its angle of inclination from the  axis, which is composed of two parts, namely, the rotation  of the cross-section induced by the pure bending and rotation  0 due to the initial imperfection.
Here, the nonlocal elastic column is treated as an elastica and the thickness-to-length ratio is assumed to be very small such that the effect of transverse shear deformation may be neglected.By assuming that the Euler-Bernoulli hypothesis holds, irrespective of small or large deformation, the strain at a distance  from the neutral axis is given by where / is the exact bending curvature.
For the allowance of size effect, Eringen's nonlocal elasticity theory is adopted.Consequently, by neglecting the nonlocal behavior in the thickness direction, the constitutive relation for a uniaxial stress state is written in the form [3,49,54] where  is the normal stress,  is the Young's modulus, and  0  is the parameter that allows for the size effect,  is an internal characteristic length (e.g., length of C-C bond, lattice spacing, and granular distance), while  0 is a constant appropriate to each material, whose magnitude is usually identified either by matching the dispersion curves of plane waves with those of atomic lattice dynamics or by calibrating it against molecular dynamic simulation results.Generally, a conservative estimate of the small-scale parameter can be set as  0  < 2.0 nm for a single-wall carbon nanotube [55].
In view of the definition of the resultant bending moment the exact moment-curvature relation can be written as where  is the flexural rigidity and  = ∫   2  is the second moment of inertia of the cross-section.
From static consideration on an arc-element  (see Figure 1(b)), we have where the horizontal and vertical internal forces () and (), for the problem considered, are given by In view of ( 5) and ( 6) and  =   −  0 , then after differentiating (4) once with respect to the arc length , the general governing equation may be written as To facilitate the manipulation of the previous system, the following dimensionless parameters are introduced: which imply that (7) can be rewritten as where the prime denotes differentiation with respect to spatial coordinate .
For the column with an initial imperfection, its initial deflection can be assumed to be [56]  0 () =  0 sin , (10) where  0 is the column's midspan initial rise.Then, from the geometrical relationships we have and correspondingly, where  0 =  0 /.The boundary condition for the postbuckling problem at hand is [57] which can also be converted into one involving the first derivative of slope   .Bearing in mind that the internal moment at any point  is  = , then from (4), one can derive through which and along with ( 12) and ( 14), the boundary condition may be given by or The use of ( 17) yields a trivial solution.Therefore, the boundary condition given by ( 16) will be used to provide a nontrivial solution.Its nondimensionalized version is as follows where   is the rotation at the end of the column.The governing differential equations ( 9) and ( 13), in conjunction with the corresponding boundary conditions (18), constitute a nonlinear two-point boundary value problem, which is different from the classical elastica theory due to the two terms on the right-hand side of (9), arising, respectively, from the size effect and initial imperfection.Owning to the strong nonlinearity, an exact or closed-form solution for this class of problem is unknown at present.Thus, the solution of this problem can only be accomplished approximately or integrated numerically.In what follows, the homotopy perturbation method and successive approximate algorithm are adopted to seek the approximate solutions.

Approximate Solution for the Differential System
Before treating the elastica problem, the zero-through thirdorder Taylor series expansions of trigonometric functions are employed as follows: so as to reduce the complexity of the boundary-value problem defined by ( 9), (13), and ( 18) and to capture the intrinsic geometrical nonlinearity at the same time.This simplification would introduce some small errors, as will be verified later by comparison with classical elastic solution.Besides, by assuming  0 ≤ 1/500 [56], ( 13) can be simplified as By substituting the expansions ( 19) and ( 20) into ( 9) and retaining terms up to O( 3  ), one obtains where L is the linear operator which implies that L(  ) =    +  2   , while the nonlinear operator defines that N(  ) =  2 [ 3   /6 −  2   (     /2 +  2  )], and the analytical function () = − 3 (1 +  2  2 ) 0 cos .Besides, the parameter  2 is defined as  2 =  2 /(1 −  2 ).

Homotopy Perturbation
Method.Among various techniques in dealing with the nonlinear differential equation (21), the homotopy perturbation method (HPM) developed by He [58,59] is one of the most effective methods.This method does not depend upon the assumption of small parameters.The main characteristic behind this approach is that, by embedding an auxiliary parameter , HPM transforms a general nonlinear problem into an infinite number of linear problems easy to solve.Its effectiveness and accuracy have been demonstrated in the analysis of various problems [60,61].
To investigate the solution of ( 21), using the homotopy technique in topology, we first construct a homotopy with an embedding parameter  ∈ [0, 1] in the form where Θ * 0 is an initial approximation; here, for the sake of simplicity, we take it as zero.The differential equation (22) satisfies the boundary conditions It is obvious that, by varying the embedding parameter  from zero to unity, (22) approaches to the original (21) from a simple equation delineated by L(Θ) − L(Θ * 0 ) = 0.According to the HPM, we assume that the solution of ( 22) can be expanded in a power series of the embedding parameter  as Furthermore, the dimensionless parameter  2 is also expanded in a series of , namely, After the substitution of the series of ( 24) and ( 25) into ( 22) and ( 23), and splitting with respect to , the following hierarchies of linear boundary value problems are obtained: where   ( = 1, 2) are defined in Appendix A.
The zero-order approximation Θ 0 () is straightforwardly determined by solving the homogeneous (26a) as which corresponds to the most important first linearized buckling mode and implies that  0 = .Introducing the last equation, along with some trigonometric identities, into (26b) results in By solving (28), we have the following expression for Θ 1 : which is accompanied by the following expression for  2 1 : Condition ( 30) is actually equivalent to the removal of a secular term as in the classical perturbation theory [40].
It should be noted that the rotation of the nonlocal elastic column in discussion can be expressed by the following base functions [30]: Therefore, similar to the first-order approximation, one can set the coefficient of cos  0  in the th-order differential equations (26a), (26b), and (26c) to zero.This provides us with the algebraic equations for the higher-order corrections of the load parameter, namely,  2  .By repeating the procedures outlined earlier, we can find sufficient accurate approximations.
If we stop at the second-order approximate solution, then   and  2 can be given, by setting  = 1, as where Ξ  ( = 1, 2) and  1 (  ;  0 , ) are defined in Appendix B.

Successive Approximate Algorithm.
The postbuckling behavior of the nonlocal elastic column is now investigated from an alternative methodology developed by Kounadis and his colleagues [62][63][64].Here, we shall refer to this alternative method as successive approximate algorithm (SAA).The convergence, uniqueness, and upper bound error estimates of solutions derived from SAA were thoroughly established in [65].
According to SAA, the reduced homogeneous linear differential equation of ( 21), namely, in conjunction with the boundary conditions (18), is first solved.Obviously, the solution  0 has the similar form as that of (26a).Then, introducing this solution into the righthand side of ( 21) yields an inhomogeneous linear differential equation ) which, associated with the conditions (18), is served as the first approximation.Then, a straightforward manipulation yields where Λ  ( 2 ,  2 0 ;   ) ( = 1, 2) are defined in Appendix C. Note that  1 is valid provided that  ̸ =  0 ( = 1, 3, . ..).By inserting (36) into the right-hand side of ( 21), one obtains Along with the boundary condition (18), the second approximate solution can be determined.By repeating this procedure, more accurate results can be reached.But in general, the first or second approximate is usually sufficient for establishing a large part of the postbuckling path [63].In view of this and considering that higher approximations require considerable computational efforts, we choose the first approximate solution (36) as the final result, which, by applying the condition  1 (0) =   , furnishes us the relation between   and  2 in the following form: Here, the coefficients , , and  are presented in Appendix C. The lowest root of last equation is given by Equation ( 39) provides us with the functional relationship of  2 versus   for various initial imperfections and smallscale parameters.
Equations ( 33) and ( 39) define the postbuckling equilibrium path.For a given rotation   at its left end, the equilibrium load  is determined by in which   =  2 / 2 is the buckling load for the same structure via local elasticity.For the perfect nonlocal elastic column, let   approach zero; then, one ends up with the critical load  cr for the onset of buckling as follows: which is an analytical solution without any approximation, and it is identical with that given in [66].Obviously, the critical load  cr is a decreasing function with respect to increasing small-scale parameter .
Once the functional relationship between   and  is known, the expected  and  coordinates of any point along the deflected neutral axis of the column can be determined by

Numerical Results and Discussion
In order to ascertain the accuracy and the range of applicability of the theoretical results developed previously, a parallel model for perfect column via local elasticity theory degenerated by setting  0 = 0 and  = 0 is first evaluated numerically against the exact elastica solution available in the literature.The comparisons of these results are presented  (a) Results by the exact theory [12]. (b) Results by the HPM truncated to the first order. (c) Results by the HPM truncated to the second order. (d) Results by the SAA.  a) Results by the exact theory [12]. (b) Results by the HPM truncated to the first order. (c) Results by the HPM truncated to the second order. (d) Results by the SAA.
in the tabular form.Tables 1 and 2 As it can be seen from the presented results, the buckling loads obtained from HPM and SAA agree well with the exact elliptical integral solution when the end rotation is less than 40 degrees, and the solution by the HPM truncated to the second order provides reliable results even for the end rotation up to 120 degrees, while for solutions provided by SAA, more iterations are needed to get accurate results.However, for the midspan deflections, the effectiveness of the aforementioned results getting both from HPM and SAA can be easily observed even when the end rotation amounts to 140 degrees.
In view of the foregoing discussions, the postbuckling behavior of the nonlocal elastic column will be identified from the results by HPM truncated to the second order.In fact, the almost identical results can be observed from using the SAA.Numerical results for perfect nonlocal elastic column are first presented in both tabular and graphical forms for various small-scale parameter .The results show that, at a specified end rotation, the size effect becomes more obvious as the postbuckling deformation increases (see Table 3 and Figure 2).To investigate the postbuckling behaviors, the stability of the column is also observed via the load-rotation curves.Figure 3 describes the size effect on the postbuckling path.It shows that the pitchfork bifurcation, composed of two symmetrical stable branches and an unstable equilibrium branch, occurs at the critical load  cr whatever the small-scale parameter values.Nevertheless, the small-scale parameter does have an appreciable effect of reducing the buckling load.As one can see from Figure 3, the deformation tends to be larger when compared to its local counterparts for the same magnitude of postbuckling load.
To illustrate the influence of the initial imperfection, several cases with or without the size effect are discussed.The bifurcation response of the imperfect column is compared in Figures 4 and 5 with that of its perfect local counterpart.As it can be seen from Figure 4, the introduction of the imperfection breaks the internal symmetry of the problem, compared with Figure 3. Buckling occurs through a saddlenode bifurcation, which makes the critical load of the column quantitatively less apparent since the critical state   is represented by the point of zero slopes on each curve.Although, for this column, the critical states other than the one for the perfect column cannot be reached under load control, it is interesting to note that the critical states of the imperfect column occur at loads higher than the critical load for its perfect counterpart, and the larger the imperfection, the greater the critical load.From Figure 4, it can also be seen that the postbuckling behavior and growth of the end rotation are altered even for a seemingly small imperfection, particularly in the neighborhood of the critical load of the perfect system, within where any slight increase of the amplitude of the imperfection would bring about greater deformation for the same load.Even so, all postbuckling paths of the imperfect system will eventually converge to the symmetrical postbuckling path of its perfect counterpart.
Unlike for the local elastic column, the postbuckling path for nonlocal elastic column additionally depends on the smallscale parameter, but the general trend of which is rather similar qualitatively as its local counterpart (see Figure 5).

Concluding Remarks
In this study, a semianalytical treatment for calculating the large elastic deformation of an initially imperfect nonlocal elastic column is presented.Herein, the column is considered to be a prismatic and inextensible one, whose constitutive equation corresponds to a differential type of Eringen's nonlocal elasticity theory.Moreover, the Euler-Bernoulli assumption is adopted.The described problem results in a complicated two-point boundary value problem with a strong nonlinearity and size effect incorporated.This problem completely precludes the use of elliptical integrals as a viable method of solution.The load-rotation relation in an explicit form, as well as the deformed curve, is obtained by the homotopy perturbation method and the successive approximate algorithm with a few iterations.Presently computed values of the postbuckling deformation and corresponding load are found to agree very well with those elastic results available in the literature.Parameter study reveals that the size effect, when the size of the column is scaled down to the nanodomains, and the initial imperfection can influence the postbuckling behavior of a nanocolumn considerably.In general, an increase in the small-scale parameter gives rise to an increase in postbuckling deformation and a decrease in the buckling load.Also, the greater the deformation becomes, the more prominent the size effect is demonstrated.Besides, the appearance of the imperfection breaks the postbuckling path from the form of an internal symmetrical pitchfork bifurcation into one of a saddle-node bifurcation.The postbuckling paths are affected primarily in the near-buckling regime; eventually, all of them will converge to its perfect counterpart.These findings will contribute to our better understanding of the special behavior of nanostructures.
From the effectiveness and accuracy of the proposed methods, we can also conclude that the presented methods can be potentially extended to a broad range of column problems under large deformations, such as the postbuckling problems of shallow arches subjected to lateral loads, problems for columns with initial imperfection having the shape of the second, or higher buckling modes [31], and problems for columns with the inextensibility assumption relaxed to an extensible one.

Figure 1 :
Figure 1: Initially imperfect nonlocal elastic column subjected to a conservative force: (a) geometry and coordinate system and (b) an infinitesimal element.

Figure 2 :Figure 3 :
Figure 2: Equilibrium configurations of a perfect nonlocal elastic column for various end rotations and small-scale parameters.

Figure 4 :Figure 5 :
Figure 4: Influence of initial imperfection on the postbuckling equilibrium paths with size effect precluded.

Table 1 :
Comparison of analytical approximations with the exact one for buckling loads.

Table 2 :
Comparison of analytical approximations with the exact one for midspan deflections.
Number inside the bracket ( ) is the % relative error computed.

Table 3 :
Midspan deflection of perfect nonlocal elastic column for various values of small-scale parameter.