Implementation of peridynamic beam and plate formulations in finite element framework

Peridynamic (PD) theory is a new continuum mechanics formulation introduced to overcome the limitations of classical continuum mechanics such as predicting crack initiation and propagation, and capturing nonlocal effects. PD theory is based on integro-differential equations and these equations are generally difficult to be solved by using analytical techniques. Therefore, numerical approximations, especially with meshless method, have been widely used. Numerical solution of three-dimensional models is usually computationally expensive and structural idealization can be utilized to reduce the computational time significantly. In this study, two of such structural idealization types are considered, namely Timoshenko beam and Mindlin plate, and their peridynamic formulations are briefly explained. Moreover, the implementation of these formulations in finite element framework is presented. To demonstrate the capability of the present approach, several case studies are considered including beam and plate bending due to transverse loading, buckling analysis and propagation of an initial crack in a plate under bending loading.


Introduction
Peridynamic (PD) theory was introduced by Silling [33] to overcome limitations of Cauchy's classical continuum mechanics (CCM) formulation for certain problems and conditions such as predicting crack initiation and propagation, and capturing nonlocal effects. PD theory is basically a new formulation of continuum mechanics. As opposed to partial differential equations of CCM, PD theory is based on integro-differential equations without having any spatial derivatives. Therefore, these equations are always valid regardless of displacement discontinuities due to cracks, etc. The derivation of peridynamic equations using Lagrange's equation is presented in Madenci and Oterkus [22]. Moreover, PD theory has an internal length scale parameter called horizon. By using horizon, it may be possible to capture nonlocal effects which can appear in various problems especially at small scales. As argued by dell'Isola et al. [5][6][7][8], the origins of PD go back to Piola's continuum formulation. Since its introduction in 2000, PD theory has been applied to many different material systems [12,18,27] and extended to other fields for heat transfer analysis, diffusion, etc. [10,16,28,29]. Queiruga and Moridis [32] performed numerical experiments on the convergence properties of state-based peridynamic laws and influence functions in two-dimensional problems.
It is difficult to solve integro-differential equations of peridynamics by using analytical techniques although analytical solutions for some simplified cases are available [24]. Instead numerical approaches are widely used mainly utilizing meshless discretization concept [34]. Meshless approaches are especially useful for analyzing large deformations. In PD theory, each material point has a domain of influence, i.e., horizon and the material point interacts with all other material points inside its horizon. Therefore, as opposed to the nearest neighbor interaction assumption of CCM, points which are far from each other can also interact in a nonlocal fashion. Moreover, in PD theory, each point has three translational degrees of freedom for three-dimensional models. To obtain the displacements of each point, three equations of motion in x-, y-and z-directions in the form of integro-differential equations should be solved. Static solution can be obtained either by using adaptive dynamic relaxation (ADR) technique [19] or solving a system of equations. For dynamic solution, after obtaining acceleration of a point from the equation of motion at particular time, implicit or explicit time integration schemes can be used to determine the new location and velocity of the point at the next time step [22].
Numerical solution of PD equations is very suitable to be done in a parallel programming environment. The solution domain can be split into multiple parts and each part can be solved by a different processor which can decrease the computational time significantly. However, even with the current high performance computing capabilities, solution of three-dimensional problems with large number of degrees of freedom and small-time step sizes can still be challenging. To overcome this problem, structural idealizations such as beam and plate formulations can be utilized. Such formulations are available in PD theory framework. For instance, Diyaroglu et al. [13] proposed an ordinary state-based Euler beam formulation to be used for slender beam structures. Taylor and Steigmann [36] developed a peridynamic plate model based on bond-based formulation by using an asymptotic analysis for thin plates. Moreover, O'Grady and Foster [25,26] developed a non-ordinary statebased peridynamic model for Euler-Bernoulli beam and Kirchhoff-Love plate formulations by disregarding the transverse shear deformations. For thick beams and plates, Diyaroglu et al. [11] introduced peridynamic Timoshenko beam and Mindlin plate formulations considering transverse shear deformations, which will also be the foundation of the current study.
Computer implementation of peridynamic formulations is usually done by using a particular programming language such as Fortran, C++. Although this approach is useful for research purposes, the applicability of these codes in a wider public may be limited since many researchers and engineers are using commercial software with user-friendly graphical user interfaces, efficient solution algorithms, etc. Therefore, in this study, the implementation of peridynamic beam and plate formulations in a commercial finite element (FE) software will be presented. Several benchmark problems are considered to demonstrate the accuracy and practicality of the current approach. Peridynamic solutions are compared against regular FE solutions [2] due to its simplicity and accuracy for problems without discontinuity [1,3,14].

Peridynamic Timoshenko beam and mindlin plate formulation
PD theory is based on the same continuum assumption as in CCM. According to this assumption, the body is composed of infinitely small volumes called material points. In the undeformed state of the body, each material point is identified by its coordinates, x (k) with (k = 1, 2, . . . , ∞), and is associated with an incremental volume, V (k) , and a mass density of ρ(x (k) ). In this section, peridynamic Timoshenko beam and Mindlin plate formulations will be briefly presented based on the study by Diyaroglu et al. [11].

Peridynamic Timoshenko beam formulation
For the material point k in a Timoshenko beam shown in Fig. 1, the peridynamic interaction forces between material points j and k arising from transverse shear deformation,f (k)( j) and bending,f (k)( j) can be defined for a linear material behavior asf in which c s and c b are the peridynamic material constants associated with the transverse shear deformation and bending of the beam, respectively. Note that the unit of the force parameters given in Eqs. (1) and (2) is actually force per volume squared since in peridynamics the equations of motion are generally written in terms of force density and the peridynamic forces are part of a volume integration. The transverse shear angle, ϕ (k)( j) , arising from the interaction between material points j and k can be defined as the average of the transverse shear angles at these material points in the form where w and φ represent the out-of-plane deflection and rotation of sections, respectively, sgn(·) function provides the sign of the associated function, and ξ ( j)(k) represents the initial distance between material points j and k. The curvature, κ (k)( j) , between the material points j and k can be expressed as Utilizing the definitions given in Eqs. (1)(2)(3)(4), the peridynamic equations of motion for a Timoshenko beam can be obtained as whereb (k) andb (k) represent body load and body moment terms, ρ, I , and A denote mass density, moment of inertia and cross-sectional area of the beam, respectively, and V ( j) is the volume of material point j.
If peridynamic interactions are limited within the horizon of a material point, then these equations can be written in integral form as and The PD material constants c s and c b can be expressed in terms of shear and Young's moduli, G and E, as where k s is the shear correction factor.

Peridynamic mindlin plate formulation
Considering the material point k as the point of interest as shown in Fig. 2, the transverse shear angle, ϕ (k)( j) , between material points j and k can be defined as the average of the shear angles at these material points in the form where φ ( j) and φ (k) represent the rotations with respect to the line of action between the material points j and k and can be defined as with θ being the orientation of the bond with respect to x-axis, φ x and φ y are rotations around y-and x-axis, respectively. As for the Timoshenko beam, the curvature, κ (k)( j) , with respect to the line of action between the material points j and k can be defined as When considering the material point j as the point of interest, the transverse shear angle and curvature for the interaction between the material points j and k become By using the same peridynamic interaction forces between material points j and k given in Eqs. (1, 2) for a Timoshenko beam and utilizing the definitions given in Eqs. (8)(9)(10)(11), the peridynamic equations of motion for a Mindlin plate can be written as or they can be expressed in continuous form as Note that the peridynamic interactions can be restricted within the horizon of material points, H . Therefore, the PD material constants c s and c b can be expressed in terms of Young's modulus, E, as where h is the thickness of the plate and Poisson's ratio is equal to ν = 1/3. The PD material parameters c b and c s are determined for a material point whose horizon is completely embedded in the material. For material points close to boundaries require correction and the details of the correction procedure are given in Diyaroglu et. al. [11].
In order to include failure in the material response, the response functions in the equations of motion for the beam and plate can be modified through a history-dependent scalar value function, where Please note that although it is possible to make connections to some of the existing approaches for fracture modeling [17,23,31,35], there are also certain differences between peridynamics and other existing techniques, especially the peridynamic bond concept and its breakage under certain conditions.

Implementation of peridynamic formulations in finite element framework
As mentioned earlier, analytical solution of peridynamic equations of motion is usually not possible and numerical approaches are widely utilized. If meshless method is used for spatial discretization, the solution domain is divided into finite number of volumes and each volume is represented by a point located at its center as shown in Fig. 3. Each point is interacting with finite number of points located inside its horizon. Peridynamic formulations can be numerically implemented using different computer languages in either serial or parallel manner. Moreover, commercial finite element softwares such as ANSYS, Abaqus can also be utilized by following the approach presented in Macek and Silling [21]. In this study, the implementation procedure of peridynamic Timoshenko beam and Mindlin plate formulations in a commercial finite element software will be presented.

Calibration process to link PD and classical FE parameters
As mentioned earlier, the unit of the peridynamic interaction force parameters between material points j and k given in Eqs. (1) and (2) is force per volume squared. In order to convert these quantities to force, these expressions should be multiplied by the volumes of the associated material points k and j as

Fig. 3 Meshless discretization of a domain and interaction between points inside the horizon of point k
In finite element framework, each peridynamic interaction can be represented using a Timoshenko beam element. Corresponding forces for a Timoshenko beam element can be expressed as in which E (k)( j) and G (k)( j) represent the Young's modulus and shear modulus of the element, respectively, and I and A are the moment of inertia and cross-sectional area. Equating the corresponding forces given in Eqs. (17) and (18) yields and Note that Young's and shear moduli expressions given in Eqs. (19a,b) are not real material property values. Instead they are serving as the calibration parameters between peridynamics and finite element framework.

Implementation in ANSYS
As explained in the previous section, a peridynamic bond can be represented using a Timoshenko beam element after calibrating the material parameters given in Eqs. (19a,b) for both Timoshenko beam and Mindlin plate formulations. In this study, ANSYS, a commercial finite element software, is utilized. ANSYS has an extensive library of finite elements with a name defined as a combination of an element type and a unique element number. A suitable element for the purpose of the current study is BEAM188. For the numerical implementation, first finite element nodes are created to represent material points. Then beam elements are created between material points to represent peridynamic interactions as shown in Fig. 4. Since in peridynamic theory a material point is interacting with other material points inside its horizon, a network of beam elements is generated to represent the numerical model as depicted in Fig. 5. If peridynamic failure criteria given in Eqs. (16a, b) are satisfied for a particular element, the element can be considered as broken utilizing EKILL command of ANSYS. Note that BEAM188 has 6 degrees of freedom. Therefore, it can both represent bending behavior and in-plane (membrane) behavior. Hence, it can capture the effect of in-plane loading on bending deformations as in the buckling analysis.

Application of boundary conditions and loading
Peridynamic equations given in Eqs. (5) and (12) are obtained using Lagrange's equations and without considering the boundaries. Although there are currently different techniques available in the literature, in this study, the displacement and rotation boundary conditions are applied by introducing a fictitious layer with a thickness equivalent to the size of horizon [22]. Moreover, the loading is applied to a single layer of material points as a body load or moment.

Numerical results
In order to demonstrate the capability of the current approach, bending and buckling analyses are performed for both peridynamic Timoshenko beam and Mindlin plate formulations. The peridynamic solutions are compared against classical finite element method solutions. All PD and FE solutions were obtained by using ANSYS software. Moreover, a plate with an initial crack under bending loading is studied to demonstrate how cracks can propagate in the current approach. The horizon size was chosen as δ = 3.015 x in all cases where x is the uniform grid spacing.

Beam bending
A cantilever beam shown in Fig. 6 with a length of L = 1 m and a cross-sectional area of A = 0.1 × 0.1 m 2 is considered. The Young's modulus and Poisson's ratio are specified as E = 200GPa and ν = 1/3, respectively. The beam is discretized into a single row of material points with a distance between each other of x = 0.01 m. As suggested by Madenci and Oterkus [22], the left edge is constrained by introducing a fictitious region with a size of horizon, δ. Therefore, a total of 103 nodes are used in the model. A transverse concentrated force P = −1000 N is applied to the right end of the beam. Figure 7 shows the deflection and rotation results obtained from both PD and classical FE models. It can be clearly seen that the solutions from PD model are in very good agreement with classical FE solution. This verifies the capability of PD Timoshenko beam formulation to capture small bending deformations and rotations.

Beam buckling
Same cantilever beam model considered in the previous case is tested for its buckling performance. When a structure is subjected to a compressive load, buckling may occur. The critical buckling load is the maximum load which a structure can resist buckling to occur.
ANSYS provides two techniques for buckling analysis which are eigenvalue analysis and geometrical nonlinear analysis. Eigenvalue analysis provides the buckling load and mode shape for each buckling mode. Note that buckling mode only demonstrates the shape of the deformation. On the other hand, geometrical nonlinear analysis can provide both the buckling load and post-buckling behavior of the structure with real deformations.

Eigenvalue buckling analysis solution
According to Euler's formula, the theoretical critical buckling load for the cantilever beam subjected to fixed-free end conditions can be calculated as The same problem is also solved by using classical finite element method and peridynamic model performing an eigenvalue analysis. Both PD and classical FE solutions of critical buckling load from eigenvalue buckling analysis yield the same value of 4.08 × 10 6 N which agree well with the theoretical value from Eq. (20). This clearly demonstrates that PD model of Timoshenko Beam has a good capability in eigenvalue analysis to predict critical buckling load.

No-linear buckling analysis solution
The performance of geometrical nonlinear buckling analysis of PD Timoshenko Beam is also studied. A total load of F x = − 5 × 10 6 N is gradually applied to the beam at the right free end using 100 load steps. A sufficiently small transverse loading of F y = − 20, 000 N(0.4% of F x ) is introduced to the free end to trigger the buckling behavior once the critical buckling load is reached. In Fig. 8, the variation of deflections as the load increases is demonstrated. The solutions of PD and classical FE agree very well with each other and both predicts that the beam becomes unstable and buckles at a load of approximately 4 × 10 6 N, which is slightly less than the eigenvalue solution of 4.08 × 10 6 N and theoretical value of 4.11 × 10 6 N. Figure 8 also shows the post-buckling behavior of the column.

Plate
In this section, the finite element implementation of PD Mindlin Plate formulation is investigated. As for the Timoshenko beam formulation, the PD Mindlin Plate formulation is tested for its bending and buckling analyses performance which include eigenvalue analysis and geometrical nonlinear analysis. PD results are compared against regular FE models created using SHELL181 element of ANSYS. Moreover, a plate with an  initial crack under bending loading is studied to demonstrate the failure prediction capability of the current formulation.

Bending
A cantilever plate with a length and width of L = W = 1 m and thickness of h = 0.1 m is considered as shown in Fig. 9. The Young's modulus and Poisson's ratio of the plate is E = 200 GPa and ν = 1/3, respectively. The model is discretized into one single row of material points along the thickness and the distance between material points is x = 0.01 m. The left edge is constrained by introducing a fictitious region with a size of horizon, δ. A transverse load is applied to a single row of material points as a body load at the right edge of the plate with an amount of 10 7 N/m 3 .
As shown in Fig. 10, the PD and classical FE solutions yield similar variation in terms of deflection and rotation for the points located along the central x-axis. This shows that FE implementation of PD can accurately capture the bending behavior of the Mindlin plate.

Buckling
The cantilever plate considered in the previous case with same dimensions and material properties is also studied for its buckling behavior. As in the Timoshenko beam case, both eigenvalue and nonlinear buckling analysis are performed in this section.

Eigenvalue buckling analysis solution
The critical buckling load of the cantilever plate is first determined by performing an eigenvalue buckling analysis. Both PD and classical FE analyses yield same critical buckling load which is approximately equal to 4.35 × 10 5 N.
Nonlinear buckling analysis solution The performance of nonlinear buckling analysis of PD Mindlin plate is studied next. A total pressure of P x = −6 × 10 8 N is gradually applied to the free edge of the plate, which is equivalent to a concentrated force of F x = −6 × 10 5 N, using 2000 load steps. To trigger the buckling behavior, a sufficiently small transverse load of F y = 2400 N (0.4% of F x ) is introduced to a single row of material points at the free edge.
In order to observe the buckling behavior of the plate, the deflection of the central point at free edge (x = L , y = 0) is recorded as the load increases. As depicted in Fig. 11, the PD and classical FE solutions are in very good agreement and both predicts that the plate buckles at a pressure loading of approximately 4.1 × 10 8 Pa which is slightly less than the eigenvalue solution of 4.35 × 10 8 Pa. Figure 11 also shows the post-buckling behavior of the cantilever plate.

Plate with an initial crack subjected to pure bending loading
As mentioned earlier, crack propagation prediction is one of the major strengths of PD theory. To demonstrate this capability, a square plate with an initial central crack aligned horizontally is considered as shown in Fig. 12. The length and width of the square plate are L = W = 1 m with a thickness of h = 0.1 m. The length of   [15], respectively. The model is discretized into one single row of material points in the thickness direction. The distance between material points is x = 0.01 m. The horizon size is chosen as δ = 3.015 x. A small increment of bending moment loading is applied through a single row of material points at the horizontal boundary regions of the plate. Under the applied uniform bending, the crack starts to propagate when the resultant body load reachesb y = ± 3.2 × 10 6 N/m. This result is consistent with the result obtained by Diyaroglu et al. [11], and as expected, the crack propagates toward the edges of the plate as the load increases as shown in Fig. 13.

Final remarks
The main purpose of this study is to present finite element implementation of peridynamic Timoshenko and Mindlin plate formulations. The advantage of this approach is that only one single row of material points along the thickness is required, which not only decreases the memory consumption by reducing the number of the nodes and elements, but also brings efficiency on processing speed. The feasibility and accuracy of the current approach is verified by considering various benchmark problems and comparing peridynamic results against classical finite element solutions in bending and buckling cases. A good agreement is obtained between peridynamic and finite element analyses results. Moreover, crack growth in a plate subjected to bending loading case is studied to demonstrate the failure prediction capability of the current approach. As a future study, impact analysis will be considered to extend the usage of the current approach. Developed framework can also be used in other applications such as bone mechanics [20]. Moreover, utilizing variational approach as presented in dell'Isola and Placidi [9] and Placidi et al. [30], the current formulation will be extended to represent the boundary conditions without utilizing fictitious regions.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.