Beam buckling analysis in peridynamic framework

Peridynamics is a non-local continuum theory which accounts for long-range internal force/moment interactions. Peridynamic equations of motion are integro-differential equations, and only few analytical solutions to these equations are available. The aim of this paper is to formulate governing equations for buckling of beams and to derive analytical solutions for critical buckling loads based on the nonlinear peridynamic beam theory. For three types of boundary conditions, explicit expressions for the buckling loads are presented. The results are compared with the classical results for buckling loads. A very good agreement between the non-local and the classical theories is observed for the case of the small horizon sizes which shows the capability of the current approach. The results show that with an increase of the horizon size the critical buckling load slightly decreases for the fixed overall stiffness of the beam.

Peridynamic formulations for elementary structures including rods and beams [3,9,23], plates [10,13,[24][25][26] and shells [4,12,17] are available. Based on the CCM, the theories of beams, plates and shells provide an efficient and robust approach for the structural analysis. Indeed, many analytical solutions for beams, plates and shells are derived in the literature providing a direct insight into the deformation and stress states. On the other hand, analytical solutions were mostly applied to validate general purpose numerical techniques and computer codes, such as the finite element method. However, within the framework of PD only few analytical solutions are available.
Analytical solutions for one-dimensional rods in both peristatic and peridynamic cases are derived in [8,15,19,21]. For statically loaded beams, truncated PD solutions are presented in [3]. The deflection function is represented with the Taylor series. Then the PD integral equation for the bending moment is solved approximately with two terms in the series expansion. This approach leads to a strain-gradient type non-local beam theory with higher derivatives of the unknown deflection function. Although several analytical solutions for PD beam equations are derived [22], explicit expressions for beam buckling loads based on PD are not presented in the literature.
The aim of this study is to derive the governing equations for buckling of beams based on the nonlinear peridynamic beam theory. For three types of boundary conditions explicit expressions for the buckling loads are presented. By use of trigonometric series representation, general solutions for the deflection function for three buckling cases are derived. Expressions for the critical buckling loads will be presented in the closed analytical form. The results will be validated by comparing the derived solutions against solutions to the classical beam theory (CBT).

Governing equations
According to CBT, the axial strain ε can be expressed as follows [6] where x is the axial coordinate, z is the transverse coordinate with origin in the cross section centroid, is the axial strain, ε is the strain of the beam centerline, u is the axial displacement, w is the deflection and κ is the curvature. Employing the Hooke's law gives the stress σ where E is the Young's modulus. The resultant force can be cast by integrating Eq. (2) over the cross section as in which A represents the cross-sectional area. For the beam subjected to compressive load ( Fig. 1), the axial resultant force must balance the external axial load, such that The strain energy per unit length of the beam according to the classical beam theory is In order to obtain the corresponding PD expression, it is necessary to introduce non-local measures for the axial strain and the curvature in Eqs (5). This can be achieved by using Taylor expansion up to second-order terms for axial and transverse displacement functions about x as follows where ξ represents the coordinate difference between material point x and its certain family member. Ignoring higher-order terms and integrating each terms with respect to ξ over the PD domain with the horizon size δ results in From Eqs (7) the nonlocal axial strainε 2 and the non-local curvatureκ are formulated as follows A class of PD theories of elastic beams can be derived assuming the correspondence of the strain energy density from the classical theory defined as a function of the PD deformation measures [26]. Applied for the Bernoulli-Euler beams, this correspondence leads to the following peridynamic strain energy W PD per unit length of the beam With Eqs (9) and (8), we obtain By setting the first variation of the potential energy to zero the following PD equations for the beam can be obtained The PD axial bond force f is expressed as follows where c is the bond stiffness and s is the bond stretch. For the homogeneous axial deformation of the beam, according to the classical beam theory the axial strain is ε = const. In this case the bond stretch is related to the axial strain by As a result from Eqs (3) and (13), the following relation between the PD bond force f and the resultant force N is obtained For the beam subjected to compressive force F, it follows With Eqs. (15) and (13), Eq. (12) takes the following form where η similarly to ξ represents the coordinate difference of a material point x and its family members. Equation (16) is the governing PD equation for the beam buckling subjected to the resultant force F.

Boundary conditions
Equation (16) is applicable for material points embedded in the PD influence domain. For material points adjacent to boundaries whose PD influence domain is defective, it is necessary to introduce fictitious material region. In what follows the width of the fictitious region is double of the PD horizon size, 2δ, as shown in Fig.  2. Let us explain several types of boundary conditions for the beam in PD framework.

Simple support
Suppose the beam is subjected to simple support at the left end. According to the classical theory the deflection and the curvature must be zero, i.e., Applying the central difference for Eq. (17) 2 yields Replacing x by ξ and associating with Eq. (17) 1 provides the PD simple support condition as It indicates that simply supported end enforces in the fictitious region an anti-symmetric relation with respect to the real material domain, as shown in Fig. 3.

Clamped support
Suppose the beam is clamped at the left end. In this case both the deflection and cross section rotation must be zero, i.e., Applying central difference with respect to the slope relation gives Equations (22) imply that clamped boundary in the fictitious region provides a symmetric relation with respect to the real material domain, as shown in Fig. 4a. Regarding roller clamped boundary (Fig. 4b), Eqs. (22) reduce to,

Free boundary
The free boundary condition imposes both the zero curvature and shear force, i.e., Applying the central difference scheme with respect to zero curvature expression gives The central difference scheme applied to zero shear force gives Coupling Eq. (25) with (26) yields It can be claimed that the relation given in Eq. (25) automatically satisfies zero shear force boundary conditions. Thus, the PD free boundary condition can be formulated as In geometrical point of view, this enforces a skew symmetric relation about point (x, z) = (L , w (L)) between fictitious region and real material vicinal to free end, as shown in Fig. 5.

Simply supported beam
Consider a simply supported beam subjected to axially compressive load as shown in Fig. 6. For the analysis, the PD buckling formulation (16) is applied with the following boundary conditions Suppose a solution to Eq. (16) satisfying the boundary conditions has the form of w (x) = a n sin nπ x L In this case one can obtain which implies the PD buckling load for simply supported beam as The PD critical buckling load is obtained for n = 1 as follows

Cantilever beam
Apparently, buckling problem of a cantilever beam is identically equivalent to that of boundary conditions of roller clamped and roller simply supported, as shown in Fig. 7. In this case Eq. (16) is applied with the following boundary conditions With consideration of boundary conditions, suppose the solution has the form of w (x) = a n cos (2n − 1) π x 2L (37) Fig. 8 Clamped -rolled clamped beam subjected to axially compressive load Substituting Eq. (37) into (16) and following the similar procedure given in the previous subsection yields the PD buckling load as The PD critical buckling load emerges when n = 1 as

Clamped -roller clamped beam
Consider a beam subjected to buckling load as shown in Fig. 8. The PD buckling formulation (16) is applied with the following boundary conditions Suppose Eq. (16) admits a solution satisfying the boundary conditions of form Plugging Eq. (41) back into (16) and obeying the procedures as in the previous cases results in the PD buckling load as in which the critical load rises at n = 1 as

Analytical validation
To verify the validity of the PD buckling load for Euler beam, PD solutions are compared against the corresponding solution in classical theory with respect to three different boundary conditions. The length of the beam is chosen as L = 1 m for each case. The results are presented in Table 1. We observe that PD solutions agree well with those of classical theory. In particular, when the PD horizon size approaches zero, PD theory converges to classical theory. With an increase of the horizon size, the PD buckling loads become slightly less than those of CBT. These results verify the accurateness of the derived PD solutions for the buckling load.

Numerical validation
To investigate whether PD buckling formulation indeed can capture the deformation undergoing buckling process, Eq. (16) is verified numerically in this section. To this end the beam is discretized by n nodes with the coordinates x ( k), k = 1, 2, . . . , n along the axis. The discretized form of Eq. (16) for the k-th node is where w ( k) is the deflection of the k-th node, w (i k ) (i = 1, 2, . . . n) is the deflection vector of the i-th node within the horizon of the node k, ξ ( j)(k) = x (k) − x ( j) and ξ (k) is the length of the node k. The geometry parameters of the beam are chosen with the length L = 1 m, thickness t = 0.05 m and width b = 0.05 m, respectively. The Young's modulus is specified as E = 200 GPa. The model is discretized into one single row of 201 material points and the horizon size is δ = 0.015 m. The axial compressive load F is applied gradually from zero till buckling phenomena emerges. Meanwhile, a vertical interfering load of q = 0.1F is applied at central point for simply supported beam and clamped-clamped beam as well as at free end for cantilever. The results are presented in Table 2. We observe that buckling phenomena indeed arises when external load approaches the critical value. These results demonstrate the validation of the PD buckling formulation in the discretized form.

Conclusions
In this study, PD buckling formulation for Euler beam was presented, which was derived by using Lagrange variational principle and Taylor expansion. The PD analytical solutions for the buckling load of beams with various boundary conditions were introduced. The validation of PD buckling formulation and analytical solutions was demonstrated by comparing the critical buckling loads against the corresponding results by the classical beam theory. A very good agreement between the non-local and the classical theories is observed for the case of the small horizon sizes which shows the capability of the current approach.
The derived analytical solutions are useful to validate the numerical methods and codes developed to solve PD equations.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
We wish to confirm that there are no known conflicts of interest associated with this publication. We confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm that the order of authors listed in the manuscript has been approved by all of us. We confirm that we have given due consideration to the protection of intellectual property associated with this work and that there are no impediments to publication, including the timing of publication, with respect to intellectual property. In so doing we confirm that we have followed the regulations of our institutions concerning intellectual property. We further confirm that any aspect of the work covered in this manuscript that has involved either experimental animals or human patients has been conducted with the ethical approval of all relevant bodies and that such approvals are acknowledged within the manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.