Peridynamic Higher-Order Beam Formulation

In this study, a novel higher-order peridynamic beam formulation is presented. The formulation is obtained by using Euler-Lagrange equations and Taylor’s expansion. To demonstrate the capability of the presented approach, several different beam configurations are considered including simply supported beam subjected to distributed loading, simply supported beam with concentrated load, clamped-clamped beam subjected to distributed loading, cantilever beam subjected to a point load at its free end and cantilever beam subjected to a moment at its free end. Transverse displacement results along the beam obtained from peridynamics and finite element method are compared with each other and very good agreement is obtained between the two approaches.


Introduction
Peridynamic (PD) theory was introduced by Silling [1] to overcome the limitations of widely used classical continuum mechanics (CCM). PD equations do not contain spatial derivatives as opposed to their usage in CCM which allows its equations to be valid even if the displacement field is not continuous due to the existence of cracks. Moreover, it has a length scale parameter called horizon which does not exist in CCM, so that it can capture nonlocal effects. PD has been applied to many different material systems including metals [2], composite materials [3] and concrete [4]. Moreover, PD is not limited to structural analysis, but can also be used to analyse heat transfer [5], diffusion [6], porous flow [7], fluid flow [8], etc. An extensive review of PD research is given in Javili et al. [9]. Original PD formulation is suitable to analyse 1-Dimensional, 2-Dimensional and 3-Dimensional models by assigning translational degrees of freedom to each material point. However, for certain shapes including beams, plates and shells, such formulations can be computational expensive. Alternatively, simplified beam, plate and shell formulations with additional rotational degrees of freedom per material point can be utilised.
There are several PD formulations available in the literature suitable for beam, plate and shell structures. Amongst these, Diyaroglu et al. [10] proposed a state-based Euler-beam formulation for slender beams. This formulation was further extended by Yang et al. [11] to analyse Kirchhoff plates. O'Grady and Foster [12] developed a non-ordinary state-based model to represent the bending behaviour of Euler-Bernoulli beam. They also introduced a non-ordinary state-based peridynamic model for Kirchhoff-Love plate [13]. In order to take into account transverse shear deformation for thick beams and plates, Diyaroglu et al. [14] presented peridynamic Timoshenko beam and Mindlin plate formulations. Taylor and Steigmann [15] developed a general plate model via bond-based peridynamics and asymptotic analysis. Chowdry et al. [16] proposed a state-based peridynamic formulation for linear elastic shells.
In this study, a novel peridynamic formulation is presented for higher-order beam theory. Higher-order beam theory allows more accurate representation in the thickness direction for thick plates with respect to Timoshenko beam formulation. Moreover, it does not require shear correction factors as in Timoshenko beam models. The formulation is based on Euler-Lagrange equation and Taylor's expansion. The capability of the proposed formulation is demonstrated by considering several benchmark problems of beams subjected to various boundary and loading conditions. Peridynamic results are compared against finite element method (FEM) results.

Classical Strain Energy Density
The displacement field of any material point in a beam can be represented in terms of the displacement field of a material point along the central axis in xz plane by using Taylor expansion as (see Fig. 1) where u(x, z) and w(x, z)are the displacement components of the material point in x-and zdirections, respectively. In this study, only flexural deformations are taken into consideration. Thus, eliminating axial deformation effects and higher order terms in Eq. (1), the components of the displacement field can be expressed as: where θ, θ * and w * are introduced as three new independent variables which are defined respectively as (see Fig. 2) In order to simplify the expressions, hereafter w(x, 0), θ(x, 0), w * (x, 0)and θ * (x, 0)will be written simply as w, θ, w * and θ * , respectively. Utilising the expressions given in Eq. (2) for the displacement components, the three well-known strain-displacement relationships of the plane stress elasticity can be written as: Note that the strain-displacement relationships given in Eq. (4a-c) are much more realistic with respect to Timoshenko beam theory where the axial normal strain, ε xx , and transverse shear strain, γ xz ,vary linearly in the thickness direction and the transverse normal strain, ε zz , is equivalent to zero. This also eliminates the need for the introduction of shear correction factor which is widely used in Timoshenko beam analysis.
Assuming the material is isotropic and plane-stress condition is applicable, the stress-strain relationships can be expressed as: where E and ν represent elastic modulus and Poisson's ratio, respectively. Substituting Eq. (4a-c) into Eq. (5) yields: The strain energy density can then be defined as: Substituting Eqs. (4) and (6) into Eq. (7) results in: The average strain energy density of a material point along the central axis of the beam can be obtained by integrating the strain energy density function given in Eq.(8) through the transverse direction and divided by the thickness, h, as:

Peridynamics
Peridynamics (PD) is a non-local theory and it is different than the classical local continuum mechanics since state of each material point is not only influenced by the material points located in its immediate vicinity but also influenced by material points which are located within a region of finite radius named as "horizon", H. The equation of motion of a material point located at x can be expressed as: where ρ(x) represents the density of the material, t represents time, u(x, t), :: u x; t ð Þ and b(x, t) are the displacement, acceleration and body load of the material point located at x. In Eq. (10), f(u ′ − u, x ′ − x) defines the peridynamic interaction (bond) force between the material point located at x and another material point inside its horizon located at x ′ (see Fig. 3).
Obtaining a closed-form solution for Eq. (10) is usually not possible. Numerical techniques including meshless approach is widely used by discretising the solution domain into small volumes and assigning a point to represent the associated volume. Therefore, the PD equations of motion for a particular material point k can be written in discrete form as: where the summation takes over the family members of the material point k, N indicates the total number of family members, and V is the volume of the material point. The PD equation of motion can be obtained by using Euler-Lagrange equation as: whereu k ð Þ is the velocity of the material point k and L is the Lagrangian which is the difference between the total kinetic energy of the system, T, and total potential energy of the system, U, i.e. L = T − U.
The total potential energy stored in the body can be obtained by summing potential energies of all material points including strain energy and energy due to external loads as: Fig. 3 Peridynamic interaction force and horizon concept [17] where i (k) is the material point inside the horizon of the material point k.
Note that unlike classical continuum mechanics, in peridynamics, the strain energy density function, W (k) , of a certain material point k depends on its displacement and the relative displacement between itself and all other material points in its horizon. The strain energy density of a material point for a higher-order beam theory based on classical continuum mechanics is given in Eq. (9). This expression can be converted into a peridynamic form by transforming the spatial derivative terms to their corresponding peridynamic form as explained in Appendix A. Therefore, by performing the following transformations, the strain energy density function given in Eq. (9) can be written in PD form as: Fig. 4 Simply supported beam subjected to distributed loading 2   6  6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  7 The generalised displacement vector, u (k) , and the body load vector, b (k) , at point k in the potential energy expression given in Eq. (13) can be written, respectively, as: The total kinetic energy of the body can be obtained by summing kinetic energies of all material points as: After substituting Eq. (2 a, b) in Eq. (17), Eq. (17) can be rewritten as: The average kinetic energy of a material point along the central axis of the beam can be obtained by integrating the kinetic energy expression given in Eq.(18) through the transverse direction and divided by the thickness, h, as: Thus, the Lagrangian of the body can be written as: Finally, Euler-Lagrange equations for higher-order beam theory can be written as:   ::

Simply Supported Beam with Distributed Load
In the first example, a simply supported beam subjected to a distributed loading of 100 N/m is considered as shown in Fig. 4. Comparison of transverse displacements obtained from PD and   Fig. 5, and a very good agreement is observed between the two approaches.

Simply Supported Beam with Concentrated Load
In the second example, the simply supported beam considered in the previous example is subjected to a concentrated load of P z = 100 N acting at the centre of the beam as shown in Fig. 6. Based on the comparison given in Fig. 7, PD results agree very well with FEM results.

Clamped-Clamped Beam with Distributed Load
In the third example, the simply supported beam subjected to distributed loading considered in the first example case is subjected clamped-clamped boundary conditions as shown in Fig. 8. The distributed load is specified as 100 N/m.
As demonstrated in Fig. 9, a very good agreement is observed between PD and FEM results in terms of transverse displacements.

Cantilever Beam Subjected to a Point Load at its Free End
In the fourth example case, a cantilever beam is considered as shown in Fig. 10. The beam is subjected to a point load of P z = 100 N at its free end. As depicted in Fig. 11, PD results agree very well with FEM results.

Cantilever Beam Subjected to a Moment at its Free End
In the final example case, the cantilever beam is subjected to a moment of 100 Nm at its free end as shown in Fig. 12. The moment loading is applied through a body load of b θ = 10 7 N/m 2 acting on a single material point at the right edge. As depicted in Fig. 13, there is a very good agreement between PD and FEM results for the transverse displacement along the beam.

Conclusion
In this study, a novel higher-order peridynamic beam formulation was presented. The formulation was obtained by using Euler-Lagrange equations and Taylor's expansion. To demonstrate the capability of the presented approach, several different beam configurations were considered including simply supported beam subjected to distributed loading, simply supported beam with concentrated load, clamped-clamped beam subjected to distributed loading, cantilever beam subjected to a point load at its free end and cantilever beam subjected to a moment at its free end. Transverse displacement results along the beam obtained from peridynamics and finite element method are compared with each other and very good agreement was obtained between the two approaches. Therefore, it can be concluded that the proposed methodology is capable of representing higher-order beam theory in peridynamic framework. To implement the clamped boundary condition, a fictitious boundary layer is created outside the actual material domain. The horizon size can be chosen asδ = 3Δxin which the discretization size is Δx. The clamped boundary condition constrains zero transverse displacement and zero rotation for the material point adjacent to the clamped end. In this study, this can be achieved by enforcing symmetrical displacement fields for w and w * and anti-symmetrical displacement fields for θ and θ * , respectively, to the material points in the fictitious region as opposed to the actual displacement field as (see Fig. 14): Simply supported boundary condition To implement the simply supported boundary condition, a fictitious layer is introduced outside the real material domain, whose size is again chosen to be equal to δ. From geometrical point of view, the simply supported boundary condition imposes zero transverse displacement and zero curvature for the material point adjacent to the constrained edge. In this study, this can be achieved by enforcing anti-symmetrical displacement fields for w and w * and symmetrical displacement fields for θ and θ * , respectively, to the material points in the fictitious region with respect to the actual displacement field as (see Fig. 15):  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/.