1 Introduction

Peridynamics (PD) is a non-local theory which accounts for long-range force/moment interactions [18]. In contrast to the classical continuum mechanics (CCM), the deformation gradient, its higher gradients or gradients of internal state variables are not required. PD equations of motion are integro-differential equations, unlike CCM where partial differential equations are derived. This allows the analysis of evolving discontinuities such as cracks in a natural way. Recent numerical studies illustrate the ability of peridynamics to capture complex fracture phenomena such as crack initiation [11, 14, 20], crack branching [7], crack kinking [5] and crack interaction with initial heterogeneities, such as holes and pores [1, 2, 16, 27].

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).

2 Peridynamic beam buckling theory

2.1 Governing equations

According to CBT, the axial strain \(\varepsilon \) can be expressed as follows [6]

$$\begin{aligned} \epsilon \left( x,z\right)= & {} \varepsilon \left( x\right) +z\kappa \left( x\right) ,\nonumber \\ {\varepsilon } \left( x\right)= & {} \frac{d{u}}{dx} +\frac{1}{2} \left( \frac{d{w}}{dx} \right) ^{2},\quad {\kappa \left( x\right) =-\frac{d^{2} {w}}{dx^{2} } }, \end{aligned}$$
(1)

where x is the axial coordinate, z is the transverse coordinate with origin in the cross section centroid, \(\epsilon \) is the axial strain, \(\varepsilon \) is the strain of the beam centerline, u is the axial displacement, w is the deflection and \(\kappa \) is the curvature. Employing the Hooke’s law gives the stress \(\sigma \)

$$\begin{aligned} \sigma \left( x,z\right) =E \epsilon =E{\varepsilon } \left( x\right) +zE\kappa \left( x\right) , \end{aligned}$$
(2)

where E is the Young’s modulus. The resultant force can be cast by integrating Eq. (2) over the cross section as

$$\begin{aligned} N \left( x\right) =\int \limits _{A}^{}\sigma dA =EA{\varepsilon } \left( x\right) \end{aligned}$$
(3)

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

$$\begin{aligned} N =EA{\varepsilon } =-F \end{aligned}$$
(4)

The strain energy per unit length of the beam according to the classical beam theory is

$$\begin{aligned} W_{\mathrm {CBT}}(\varepsilon ^2,\kappa ^2) {=} \frac{1}{2} EA\varepsilon ^2+ \frac{1}{2} EI\kappa ^2 \end{aligned}$$
(5)

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

$$\begin{aligned} {u}\left( x+\xi \right) -{u}\left( x\right) =\frac{\partial {u}}{\partial x} \xi +\frac{1}{2} \frac{\partial ^{2} {u}}{\partial x^{2} } \xi ^{2}, \quad {w}\left( x+\xi \right) -{w}\left( x\right) =\frac{\partial {w}}{\partial x} \xi +\frac{1}{2} \frac{\partial ^{2} {w}}{\partial x^{2} } \xi , ^{2} \end{aligned}$$
(6)

where \(\xi \) represents the coordinate difference between material point x and its certain family member. Ignoring higher-order terms and integrating each terms with respect to \(\xi \) over the PD domain with the horizon size \(\delta \) results in

$$\begin{aligned} \begin{array}{l} \displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left\{ u\left( x+\xi \right) -u\left( x\right) +\frac{1}{2} \frac{\left[ w\left( x+\xi \right) -w\left( x\right) \right] ^{2} }{\xi } \right\} ^{2}d\xi =\left[ \frac{du}{dx} +\frac{1}{2} \left( \frac{dw}{dx} \right) ^{2} \right] ^{2} \int _{-\delta }^{\delta }\left| \xi \right| d\xi , \\ \displaystyle \int \limits _{-\delta }^{\delta }\frac{{w}\left( x+\xi \right) -{w}\left( x\right) }{\xi ^{2} } d\xi =\frac{\partial {w}}{\partial x} \int \limits _{-\delta }^{\delta }\frac{1}{\xi } d\xi +\frac{1}{2} \frac{\partial ^{2} {w}}{\partial x^{2} } \int \limits _{-\delta }^{\delta }d\xi \end{array} \end{aligned}$$
(7)

From Eqs (7) the nonlocal axial strain \(\bar{\varepsilon }^2\) and the non-local curvature \(\bar{\kappa }\) are formulated as follows

$$\begin{aligned} \begin{array}{l} \displaystyle \bar{\varepsilon }^2=\frac{1}{\delta ^{2} } \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left\{ u\left( x+\xi \right) -u\left( x\right) +\frac{1}{2} \frac{\left[ w\left( x+\xi \right) -w\left( x\right) \right] ^{2} }{\xi } \right\} ^{2} d\xi , \\ \displaystyle \bar{\kappa }=\frac{1}{\delta } \int \limits _{-\delta }^{\delta }\frac{w\left( x+\xi \right) -w\left( x\right) }{\xi ^{2} } d\xi \end{array} \end{aligned}$$
(8)

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_{\mathrm {PD}}\) per unit length of the beam

$$\begin{aligned} W_{\mathrm {PD}}=W_{\mathrm {CBT}}(\bar{\varepsilon }^2,\bar{\kappa }^2) \end{aligned}$$
(9)

With Eqs (9) and (8), we obtain

$$\begin{aligned} W_{\mathrm {PD}}= \frac{EA}{2\delta ^{2} } \int \limits _{-\delta }^{\delta }\frac{1}{|\xi |} \left\{ u(x+\xi )-u(x)+\frac{1}{2} \frac{\left[ w(x+\xi )-w(x)\right] ^{2} }{\xi } \right\} ^{2} d\xi +\frac{EI}{2\delta ^{2} } \left( \int \limits _{-\delta }^{\delta }\frac{w(x+\xi )-w(x)}{\xi ^{2} } d\xi \right) ^{2}\nonumber \\ \end{aligned}$$
(10)

By setting the first variation of the potential energy to zero

$$\begin{aligned} \delta \Pi = 0, \quad \Pi =\int \limits _0^L W_{\mathrm {PD}}(x)\mathrm {d} x- \int \limits _0^L p(x){w}(x)\mathrm {d} x, \end{aligned}$$
(11)

the following PD equations for the beam can be obtained

$$\begin{aligned} \begin{array}{l} \displaystyle \frac{2EA}{\delta ^{2} } \int \limits _{-\delta }^{\delta }\frac{1}{|\xi |} \left\{ u(x+\xi )-u(x)+\frac{1}{2} \frac{\left[ w(x+\xi )-w(x)\right] ^{2} }{\xi } \right\} d\xi =0,\\ \displaystyle \frac{2EA}{\delta ^{2} } \int \limits _{-\delta }^{\delta }\frac{1}{|\xi |} \frac{w(x+\xi )-w(x)}{\xi } \left\{ u(x+\xi )-u(x)+\frac{1}{2} \frac{\left[ w(x+\xi )-w(x)\right] ^{2} }{\xi } \right\} d\xi \\ \displaystyle \qquad +\frac{EI}{\delta ^{2} } \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left[ \int _{-\delta }^{\delta }\frac{w(x+\eta )-w(x)}{\eta ^{2} } d\eta -\int _{-\delta }^{\delta }\frac{w(x+\xi +\eta )-w(x+\xi )}{\eta ^{2} } d\eta \right] d\xi +p(x)=0 \end{array} \end{aligned}$$
(12)

The PD axial bond force f is expressed as follows

$$\begin{aligned} f=c s, \quad c=\frac{2EA}{\delta ^{2} }, \quad s=\frac{1}{|\xi |} \left\{ u(x+\xi )-u(x)+\frac{1}{2} \frac{\left[ w(x+\xi )-w(x)\right] ^{2} }{\xi } \right\} , \end{aligned}$$
(13)

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 \(\varepsilon =\mathrm {const}\). In this case the bond stretch is related to the axial strain by

$$\begin{aligned} s\, \mathrm {sgn}\xi =\varepsilon \end{aligned}$$

As a result from Eqs (3) and (13), the following relation between the PD bond force f and the resultant force N is obtained

$$\begin{aligned} f\mathrm {sgn}\xi =\frac{2N}{\delta ^2} \end{aligned}$$
(14)

For the beam subjected to compressive force F, it follows

$$\begin{aligned} f\mathrm {sgn}\xi =-\frac{2F}{\delta ^2} \end{aligned}$$
(15)

With Eqs. (15) and (13), Eq. (12) takes the following form

$$\begin{aligned} \begin{array}{l} \displaystyle -\frac{2F}{\delta ^{2} } \int \limits _{-\delta }^{\delta }\frac{{w}\left( x+\xi \right) -{w}\left( x\right) }{\left| \xi \right| } d\xi \\ \displaystyle \quad +\frac{EI}{\delta ^{2} } \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left[ \int \limits _{-\delta }^{\delta }\frac{{w}\left( x+\eta \right) -{w}\left( x\right) }{\eta ^{2} } d\eta -\int \limits _{-\delta }^{\delta }\frac{{w}\left( x+\xi +\eta \right) -{w}\left( x+\xi \right) }{\eta ^{2} } d\eta \right] d\xi +p\left( x\right) =0, \end{array} \end{aligned}$$
(16)

where \(\eta \) similarly to \(\xi \) 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.

Fig. 1
figure 1

Beam under transverse and compressive loading

Full size image

2.2 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\delta \), as shown in Fig. 2. Let us explain several types of boundary conditions for the beam in PD framework.

Fig. 2
figure 2

Real and fictitious domains for the beam

Full size image

2.2.1 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.,

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle {w\left( 0\right) =0}, \\ \displaystyle {\left. \frac{\mathrm {d}^{2} w}{\mathrm {d}x^{2} } \right| _{x=0} =0} \end{array}\right. \end{aligned}$$
(17)

Applying the central difference for Eq. (17)\(_2\) yields

$$\begin{aligned} \left. \frac{\mathrm {d}^{2} w}{\mathrm {d}x^{2} } \right| _{x=0} \approx \frac{w\left( -\Delta x\right) -2w\left( 0\right) +w\left( \Delta x\right) }{\Delta x^{2} } =0 \end{aligned}$$
(18)

Replacing \(\Delta x\) by \(\xi \) and associating with Eq. (17)\(_1\) provides the PD simple support condition as

$$\begin{aligned} {w\left( -\xi \right) =-w\left( \xi \right) }, \quad \text{ for } \quad {0\le \xi \le 2\delta } \end{aligned}$$
(19)

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.

Fig. 3
figure 3

Simply supported beam with a fictitious region

Full size image

2.2.2 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.,

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle {w\left( 0\right) =0}, \\ \displaystyle {\left. \frac{\mathrm {d}w}{\mathrm {d}x} \right| _{x=0} =0} \end{array}\right. \end{aligned}$$
(20)

Applying central difference with respect to the slope relation gives

$$\begin{aligned} \left. \frac{\mathrm {d}w}{\mathrm {d}x} \right| _{x=0} \approx \frac{w\left( \Delta x\right) -w\left( -\Delta x\right) }{2\Delta x} =0 \end{aligned}$$
(21)

Replacing \(\Delta x\) by PD notation \(\xi \) and associating with Eq. (20)\(_1\) gives the PD clamped boundary condition as

$$\begin{aligned} \begin{array}{cc} {\left\{ \begin{array}{l} \displaystyle {w\left( 0\right) =0}, \\ \displaystyle {w\left( -\xi \right) =w\left( \xi \right) } \end{array}\right. }&{\qquad \text{ for }\quad 0\le \xi \le 2\delta } \end{array} \end{aligned}$$
(22)

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.

Fig. 4
figure 4

Clamped beam with a fictitious region

Full size image

Regarding roller clamped boundary (Fig. 4b), Eqs. (22) reduce to,

$$\begin{aligned} \begin{array}{cc} {w\left( -\xi \right) =w\left( \xi \right) }&{\qquad \text{ for }\quad 0\le \xi \le 2\delta } \end{array} \end{aligned}$$
(23)

2.2.3 Free boundary

The free boundary condition imposes both the zero curvature and shear force, i.e.,

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle {\left. \frac{\mathrm {d}^{2} w}{\mathrm {d}w^{2} } \right| _{x=L} =0} \\ \displaystyle {\left. \frac{\mathrm {d}^{3} w}{\mathrm {d}w^{3} } \right| _{x=L} =0} \end{array}\right. \end{aligned}$$
(24)

Applying the central difference scheme with respect to zero curvature expression gives

$$\begin{aligned} {\left. \frac{\mathrm {d}^{2} w}{\mathrm {d}w^{2} } \right| _{x=L} \approx \frac{w\left( L-\Delta x\right) -2w\left( L\right) +w\left( L+\Delta x\right) }{\Delta x^{2} } } {=} {0} \quad {\Rightarrow \quad w\left( L+\xi \right) } {=} {2w\left( L\right) -w\left( L-\xi \right) } \end{aligned}$$
(25)

The central difference scheme applied to zero shear force gives

$$\begin{aligned} \begin{array}{l} \displaystyle {\left. \frac{\mathrm {d}^{3} w}{\mathrm {d}w^{3} } \right| _{x=L} \approx \frac{-w\left( L-2\Delta x\right) +2w\left( L-\Delta x\right) -2w\left( L+\Delta x\right) +w\left( L+2\Delta x\right) }{2\Delta x^{3} } } {=} {0} \\ \displaystyle {\Rightarrow -w\left( L-2\xi \right) +2w\left( L-\xi \right) -2w\left( L+\xi \right) +w\left( L+2\xi \right) } {=} {0} \end{array} \end{aligned}$$
(26)

Coupling Eq. (25) with (26) yields

$$\begin{aligned} \begin{array}{l} \displaystyle {-w\left( L-2\xi \right) +2w\left( L-\xi \right) -2w\left( L+\xi \right) +w\left( L+2\xi \right) }\\ \displaystyle \quad {=} {-2w\left( L-2\xi \right) +4w\left( L-\xi \right) -2w\left( L\right) } {\approx -2\left. \xi ^{2} \frac{\mathrm {d}^{2} w}{\mathrm {d}x^{2} } \right| _{x=L} } \end{array} \end{aligned}$$
(27)

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

$$\begin{aligned} \begin{array}{cc} {w\left( L+\xi \right) =2w\left( L\right) -w\left( L-\xi \right) }&{\qquad \text{ for }\quad 0\le \xi \le 2\delta } \end{array} \end{aligned}$$
(28)

In geometrical point of view, this enforces a skew symmetric relation about point \(\left( x,z\right) =\left( L,w\left( L\right) \right) \) between fictitious region and real material vicinal to free end, as shown in Fig. 5.

Fig. 5
figure 5

A beam with a free boundary

Full size image

3 Analytical solutions

3.1 Simply supported beam

Fig. 6
figure 6

Simply supported beam subjected to axially compressive load

Full size image

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

$$\begin{aligned} \displaystyle \left\{ \begin{array}{ll} {{w}\left( -\xi \right) =-{w}\left( \xi \right) }, &{} {0\le \xi \le 2\delta } \\ \displaystyle {{w}\left( L+\xi \right) =-{w}\left( L-\xi \right) }, &{} {0\le \xi \le 2\delta } \end{array}\right. \end{aligned}$$
(29)

Suppose a solution to Eq. (16) satisfying the boundary conditions has the form of

$$\begin{aligned} {w}\left( x\right) =a_{n} \sin \left( \frac{n\pi x}{L} \right) \end{aligned}$$
(30)

In this case one can obtain

$$\begin{aligned} {\int \limits _{-\delta }^{\delta }\frac{{w}\left( x+\xi \right) -{w}\left( x\right) }{\left| \xi \right| } d\xi } {=} {a_{n} \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left[ \sin \frac{n\pi \left( x+\xi \right) }{L} -\sin \frac{n\pi x}{L} \right] d\xi } \end{aligned}$$
(31)

Likewise, we have

$$\begin{aligned} \begin{array}{l} \displaystyle \int \limits _{-\delta }^{\delta }\frac{{w}\left( x+\eta \right) -{w}\left( x\right) }{\eta ^{2} } d\eta =a_{n} \sin \frac{n\pi x}{L} \int \limits _{-\delta }^{\delta }\frac{1}{\eta ^{2} } \left( \cos \frac{n\pi \eta }{L} -1\right) d\eta , \\ \displaystyle \int \limits _{-\delta }^{\delta }\frac{{w}\left( x+\xi +\eta \right) -{w}\left( x+\xi \right) }{\eta ^{2} } d\eta =a_{n} \sin \frac{n\pi \left( x+\xi \right) }{L} \int \limits _{-\delta }^{\delta }\frac{1}{\eta ^{2} } \left( \cos \frac{n\pi \eta }{L} -1\right) d\eta \end{array} \end{aligned}$$
(32)

Coupling Eqs (31) and (32) with (16) results in

$$\begin{aligned} \begin{array}{l} \displaystyle {-\frac{2F}{EI} a_{n} \sin \frac{n\pi x}{L} \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left( \cos \frac{n\pi \xi }{L} -1\right) d\xi +} \\ \displaystyle {\int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left[ a_{n} \sin \frac{n\pi x}{L} \int \limits _{-\delta }^{\delta }\frac{1}{\eta ^{2} } \left( \cos \frac{n\pi \eta }{L} -1\right) d\eta -a_{n} \sin \frac{n\pi \left( x+\xi \right) }{L} \int \limits _{-\delta }^{\delta }\frac{1}{\eta ^{2} } \left( \cos \frac{n\pi \eta }{L} -1\right) d\eta \right] d\xi =0} \end{array} \end{aligned}$$
(33)

which implies the PD buckling load for simply supported beam as

$$\begin{aligned} F=\frac{EI}{2} \frac{\displaystyle \left[ \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left( 1-\cos \frac{n\pi \xi }{L} \right) d\xi \right] ^{2} }{\displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left( 1-\cos \frac{n\pi \xi }{L} \right) d\xi } \end{aligned}$$
(34)

The PD critical buckling load is obtained for \(n=1\) as follows

$$\begin{aligned} F_{\mathrm {cr}} =\frac{EI}{2} \frac{\displaystyle \left[ \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left( 1-\cos \frac{\pi \xi }{L} \right) d\xi \right] ^{2} }{\displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left( 1-\cos \frac{\pi \xi }{L} \right) d\xi } \end{aligned}$$
(35)

3.2 Cantilever beam

Fig. 7
figure 7

Cantilever beam subjected to axially compressive load

Full size image

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

$$\begin{aligned} \left\{ \begin{array}{cc} \displaystyle {{w}\left( -\xi \right) ={w}\left( \xi \right) }, &{} {0\le \xi \le 2\delta } \\ {{w}\left( L+\xi \right) =-{w}\left( L-\xi \right) } &{} {0\le \xi \le 2\delta } \end{array}\right. \end{aligned}$$
(36)

With consideration of boundary conditions, suppose the solution has the form of

$$\begin{aligned} {w}\left( x\right) =a_{n} \cos \frac{\left( 2n-1\right) \pi x}{2L} \end{aligned}$$
(37)

Substituting Eq. (37) into (16) and following the similar procedure given in the previous subsection yields the PD buckling load as

$$\begin{aligned} F=\frac{EI}{2} \frac{\displaystyle \left\{ \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left[ \cos \frac{\displaystyle \left( 2n-1\right) \pi \xi }{2L} -1\right] d\xi \right\} ^{2} }{\displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left[ \cos \frac{\left( 2n-1\right) \pi \xi }{2L} -1\right] d\xi } \end{aligned}$$
(38)

The PD critical buckling load emerges when \(n=1\) as

$$\begin{aligned} F_{\mathrm {cr}} =\frac{EI}{2} \frac{\displaystyle \left\{ \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left[ \cos \frac{\pi \xi }{2L} -1\right] d\xi \right\} ^{2} }{\displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left[ \cos \frac{\pi \xi }{2L} -1\right] d\xi } \end{aligned}$$
(39)

3.3 Clamped – roller clamped beam

Fig. 8
figure 8

Clamped – rolled clamped beam subjected to axially compressive load

Full size image

Consider a beam subjected to buckling load as shown in Fig. 8. The PD buckling formulation (16) is applied with the following boundary conditions

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle {w\left( 0\right) =0} \\ {w\left( L\right) =0} \\ {\begin{array}{ll} {w\left( -\xi \right) =w\left( \xi \right) ,} &{} {0\le \xi \le 2\delta } \\ {w\left( L+\xi \right) =-w\left( L-\xi \right) ,} &{} {0\le \xi \le 2\delta } \end{array}} \end{array}\right. \end{aligned}$$
(40)

Suppose Eq. (16) admits a solution satisfying the boundary conditions of form

$$\begin{aligned} {w}\left( x\right) =a_{n} \left( \cos \frac{2n\pi }{L} -1\right) \end{aligned}$$
(41)

Plugging Eq. (41) back into (16) and obeying the procedures as in the previous cases results in the PD buckling load as

$$\begin{aligned} F=\frac{EI}{2} \frac{\left\{ \displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left[ \cos \frac{2n\pi \xi }{L} -1\right] d\xi \right\} ^{2} }{\displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left[ \cos \frac{2n\pi \xi }{L} -1\right] d\xi } \end{aligned}$$
(42)

in which the critical load rises at \(n=1\) as

$$\begin{aligned} F_{\mathrm {cr}} =\frac{EI}{2} \frac{\displaystyle \left\{ \int \limits _{-\delta }^{\delta }\frac{1}{\xi ^{2} } \left[ \cos \frac{2\pi \xi }{L} -1\right] d\xi \right\} ^{2} }{\displaystyle \int \limits _{-\delta }^{\delta }\frac{1}{\left| \xi \right| } \left[ \cos \frac{2\pi \xi }{L} -1\right] d\xi } \end{aligned}$$
(43)

4 Validation

4.1 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.

Table 1 Critical buckling loads for several cases
Full size table

4.2 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,\ldots ,n\) along the axis. The discretized form of Eq. (16) for the k-th node is

$$\begin{aligned} \begin{array}{l} \displaystyle -\frac{2F}{\delta ^{2} } \sum _{j}^{}\frac{{w}_{\left( j\right) } -{w}_{\left( k\right) } }{\left| \xi _{\left( j\right) \left( k\right) } \right| } \Delta \xi _{\left( j\right) }\\ \quad \displaystyle +\frac{EI}{\delta ^{2} } \sum _{j}^{}\frac{1}{\xi _{\left( j\right) \left( k\right) }^{2} } \left( \sum _{i^{k} }^{}\frac{{w}_{\left( i^{k} \right) } -{w}_{\left( k\right) } }{\xi _{\left( i^{k} \right) \left( k\right) }^{2} } \Delta \xi _{\left( i^{k} \right) } -\sum _{i^{j} }^{}\frac{{w}_{\left( i^{j} \right) } -{w}_{\left( j\right) } }{\xi _{\left( i^{j} \right) \left( j\right) }^{2} } \Delta \xi _{\left( i^{j} \right) } \right) \Delta \xi _{\left( j\right) } +p_{\left( k\right) } =0, \end{array} \end{aligned}$$
(44)

where \(w_(k)\) is the deflection of the k-th node, \({w}_{\left( i^{k} \right) }\) (\(i=1, 2, \ldots n\)) is the deflection vector of the i-th node within the horizon of the node k, \(\xi _{\left( j\right) \left( k\right) }=x_{(k)}-x_{(j)}\) and \(\Delta \xi _{\left( k\right) }\) 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 \(\delta =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.

Table 2 Deflection versus applied force several cases
Full size table

5 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.