A State-Based Peridynamic Formulation for Functionally Graded
Euler-Bernoulli Beams

In this study, a new state-based peridynamic formulation is developed for functionally graded Euler-Bernoulli beams. The equation of motion is developed by using Lagrange’s equation and Taylor series. Both axial and transverse displacements are taken into account as degrees of freedom. Four different boundary conditions are considered including pinned support-roller support, pinned support-pinned support, clamped-clamped and clamped-free. Peridynamic results are compared against finite element analysis results for transverse and axial deformations and a very good agreement is observed for all different types of boundary conditions.


Introduction
Functionally Graded Materials (FGMs) are within the class of composite materials which have material properties continuously varying from one surface to another surface. This eliminates the stress concentrations that occur at the interface of neighbouring plies in laminated composites. Therefore, failure due to delamination can be avoided. FGMs are used in different fields including mechanical, aerospace, nuclear and civil engineering [1].
Due to the increase in usage of FGMs, numerous beam formulations for FGMs are available in the literature. Amongst these, Li [2] developed a unified approach to analyze the static and dynamic behaviours of functionally graded beams by including rotary inertia and shear deformation. Kadoli et al. [3] performed static analysis of functionally graded beams using higher order shear deformation theory. Kahrobaiyan et al. [4] proposed a size-dependent functionally graded Euler-Bernoulli beam model based on strain gradient theory. Thai et al. [1] performed bending and free vibration analysis of functionally graded beams using various higher-order shear deformation beam theories. Lu et al. [5] presented elasticity solutions for bending and thermal deformations of functionally graded beams with various end conditions using the state space-based differential quadrature method. Su et al. [6] developed the dynamic stiffness method to investigate the free vibration behaviour of functionally graded beams. Simsek et al. [7] performed linear dynamic analysis of an axially functionally graded beam with simply supported edges and subjected to a moving harmonic load based on Euler-Bernoulli beam theory. Li et al. [8] investigated relationships between buckling loads of functionally graded Timoshenko and homogenous Euler-Bernoulli beams. Sankar [9] presented an elasticity solution for a functionally graded Euler-Bernoulli beam subjected to transverse loads. Birsan et al. [10] performed deformation analysis of functionally graded beams by using direct approach which is based on the deformable curve model with a triad of rotating directors attached to each point. Nguyen et al. [11] investigated static and free vibration behaviour of axially loaded functionally graded beams based on the first-order shear deformation theory. Filippi et al. [12] used 1D Carrera Unified Formulation (CUF) to perform static analysis of functionally graded beams.
In this study, an alternative formulation, peridynamics [13][14][15][16][17][18][19][20][21][22], is utilized to analyze functionally graded Euler-Bernoulli beams. Peridynamics is a non-local continuum mechanics formulation and it is very suitable for failure analysis of structures due to its mathematical structure. Moreover, it has a length scale parameter called "horizon" which provides a capability to represent non-classical deformation behaviour which are usually seen in objects at nano-scale. The nonlocal operator introduced by Rabczuk et al. [23] and Ren et al. [24] can significantly simplify the derivation of the strong form and weak form of peridynamics. The novel peridynamic formulation specifically developed for functionally graded Euler-Bernoulli beams is based on state-based peridynamics and the equation of motion is obtained utilising Lagrange's equation. Several different numerical cases are considered to demonstrate the capability of the current approach for functionally graded Euler-Bernoulli beams subjected to different types of boundary conditions.

Euler-Bernoulli Beam Theory
A complete and adequate set of equation for linear theory of thin beams was developed by Euler and Bernoulli, which is also known as Euler-Bernoulli beam theory. According to Euler-Bernoulli beam theory, a transverse normal to the central axis of the beam in the undeformed state remains straight, normal to the central axis and its length doesn't change during deformation. Based on the assumptions of the Euler-Bernoulli beam theory, the displacement field of any material point can be represented in terms of the displacement field of the material points along the central axis in xz plane as where u x; 0 ð Þ and w x; 0 ð Þ denote the displacement components of the material point along the central axis in x-and z-directions, respectively. Thus, the strain-displacement relationship can be written as in which h ¼ @u 0 @z represents the rotation angle of the material point along the central axis, which can be expressed in terms of w 0 by casting from Eq. (2b) as Therefore, eliminating h, the deformation in an Euler-Bernoulli beam can be easily represented in terms of only two central axis displacement components u 0 , and w 0 as where u 0 ¼ u x; 0 ð Þ and w 0 ¼ w x; 0 ð Þ represent the displacement functions of the central axis of the beam, respectively.
According to the Hooke's Law, the stress function can be written as 3 Peridynamic Theory Peridynamic (PD) theory is a reformulation of continuum mechanics by using integro-differential equations instead of partial differential equations. The formulation does not contain any spatial derivatives which prevents issues due to discontinuities especially in the displacement field as a result of cracks. Material points in PD formulation can interact with each other in a non-local manner within a domain of influence called "horizon". The governing equations of peridynamics can be obtained by using Lagrange equations [25] and can be written as where q is the density, u is the displacement of the material point located at x, and € u is the acceleration of the same material point. f represents the interaction (bond) force between two interacting material points located at x and x 0 . b is the external body force acting on the material point located at x. Note that the material points represented with x 0 are the material points within the horizon, H of the material point located at x. To obtain a closed-form solution for Eq. (9) is usually not possible and numerical methods are generally utilized. Therefore, the discrete form of Eq. (9) for a particular material point k can be written as where N is the number of material points inside the horizon of the material point k. The symbol j represents the material points inside the horizon of the material point k and the bond force between the material points k and j can be expressed as where t k ð Þ j ð Þ is the force that the material point j exerting on k and t j ð Þ k ð Þ is the force that the material point k exerting on j as depicted in Fig. 1.
As mentioned earlier, the PD equations of motion can be obtained by using Lagrange's equation which can be written for the material point k as where the Lagrangian term, L can be defined as the difference between the total kinetic energy, T and potential energy, U of the system. The total potential energy of the system is the difference between the total strain energy of the system and the total energy due to external forces which can be expressed as where W k ð Þ is the strain energy density of the material point k and i k i ¼ 1; 2; 3; Á Á Á ð Þcorresponds to the material points inside the horizon of the material point k. On the other hand, the total kinetic energy of the system can be calculated as where _ u k ð Þ is the velocity of the material point k. Utilising Eqs. (12) and (13), the Lagrangian term can be written as Hence, by utilizing the Lagrange's equation given in Eq. (12) and the Lagrangian given in Eq. (15), PD equations of motion can be obtained as where the material point i k k is a material point inside the horizon of the material point i k . Since the material point k is inside the horizon of the material point the material point i k , by using the following definition PD equations of motion can be simplified as which is equivalent to the PD equations of motion given in Eq. (9). Hence, the interaction force expressions can be written as a function of strain energy density functions as For a functionally graded beam, the components of the PD interaction forces in x-(axial) and z-(transverse) directions can be written as and where u and w are the components of the displacement vector in x-and z-directions, respectively.
In order to express the strain energy density function for the material point k in non-local form, it is necessary to transform all differential terms to their equivalent form of integration. To achieve this, first, Taylor expansion of the axial displacement function u x ð Þ can be written while ignoring higher order terms as Squaring both sides of Eq. (21) and dividing all terms by the absolute value of n results in Considering x as a constant point and integrating both sides of Eq. (22) over a symmetric domain, Next, dividing both sides of Eq. (21) by n results in Considering x as a constant point and integrating both sides of Eq. (24) over a symmetric domain, Similarly, the transverse displacement function w x ð Þ can be expanded using Taylor series while ignoring higher order terms as Multiplying each term in Eq. (26) by 1 n 2 and integrating over a symmetric domain, Àd; d ½ , yields Note that the transverse displacement function w x ð Þ is related to the flexural deformation of the beam, thus, multiplying Eq. (26) by 1 n 2 is for the sake of ensuring the dimension of integrand of Eq. (27) in accordance to the curvature, i.e., "1/length".
Strain energy density expression can then be written in its non-local form by substituting Eqs. (23), (24), and (27) into Eq. (8) as It can also be written in discretized form for the material point k as where A represents the cross-sectional area of the beam. Similarly, the PD strain energy density function for the material point j can be written by changing index as Substituting Eqs. (29a) and (29b) into Eqs. (20a) and (20b) yields the PD force densities as The derivation procedure of Eqs. (31) and (32) is given in Appendix C. The final peridynamic equation of motion for a functionally graded Euler-Bernoulli beam can be obtained by inserting Eqs. (30) and (31) into Eq. (10) as

Boundary Conditions
As mentioned earlier, PD equations of motion are in the form of integro-differential equations which are different than the partial differential equations of Classical Continuum Mechanics (CCM). Application of boundary conditions in PD is also different than the ones in CCM and they are applied over a volume. In this section, two common boundary condition types for beams are considered and the procedure of application of such boundary conditions in PD framework is explained in detail.

Clamped Boundary Condition
According to Euler-Bernoulli Beam theory, the clamped boundary condition can be represented by constraining the transverse displacement, w and the rotation, h at the boundary which can be defined as In PD framework, such boundary condition can be achieved by introducing a fictitious domain outside the boundary with a length equal to two times of the horizon size, d so that all material points inside the actual solution domain have a complete horizon. In this study, the horizon size is chosen as d ¼ 3Dx where Dx is the distance between material points. Therefore, there are six additional material points inside the fictitious boundary region.
As shown in Fig. 2, clamped boundary conditions given in Eq. (33) can be expressed in PD form by imposing mirror images of the transverse displacements for the six material points in the actual and fictitious domains next to the boundary as . . . ; 6 (34a) and Figure 2: Application of clamped boundary condition In addition, the axial displacement should also be constrained at the clamped boundary as which can be expressed in PD framework by imposing anti-mirror values of the axial displacements of the six material points in the actual and fictitious domains as . . . ; 6 (36)

Simply Supported Boundary Condition
According to Euler-Bernoulli Beam theory, the simply supported boundary condition can be represented by constraining transverse displacement, w and curvature at the boundary which can be defined as As in the clamped boundary condition, simply supported boundary condition in PD framework can be achieved by introducing a fictitious domain outside the boundary with a length equal to two times of the horizon size, d. Moreover, as shown in Figs. 3a and 3b, the constraint conditions given in Eq. (37) can be satisfied by imposing anti-symmetrical transverse displacement fields to the six material points in the actual and fictitious domains next to the boundary which can be defined as . . . ; 6 (38) However, for axial deformation, the application of simply supported boundary condition is different depending on the pinned support and roller support conditions. For pinned support condition, the axial deformations can be constrained by imposing anti-symmetrical displacement fields to the six material points in the actual and fictitious material domains next to the boundary, as shown in Fig. 3a, which can be defined as On the other hand, the implementation of roller support boundary condition requires symmetrical displacement field adjacent to the boundary as (see Fig. 3b) . . . ; 6 (40)

Numerical Results
To verify the validity of the PD formulation for functionally graded Euler-Bernoulli beams, the PD solutions are compared with the corresponding finite element (FE) analysis results. In this study, the functionally graded material property is chosen as Young's Modulus, E z ð Þ, and it is assumed to vary linearly through the thickness as where E t and E b denote the Young's modulus of the top and bottom surfaces of the beam, and h donates the total thickness of the beam as shown in Fig. 4. The procedure to calculate surface correction factors is presented in Appendix A. All numerical examples considered in this section are for statics analysis and the numerical solution procedure is given in Appendix B. For all numerical examples, the horizon size is chosen as d ¼ 3Dx where Dx is the distance between material points.

Functionally Graded Beam Subjected to Pinned Support-Roller Support Boundary Condition
A simply supported functionally graded beam with length, width and thickness of L ¼ 1 m, b ¼ 0:01 m and h ¼ 0:05 m, is considered as shown in Fig. 5. The beam is constrained by pinned support and roller support at left and right ends, respectively. The Young's modulus of the top and bottom surfaces are E t ¼ 200 GPa and E b ¼ 100 GPa, respectively. The model is discretized into one single row of material points along the thickness and the distance between material points is Dx ¼ 1=101 m. A fictitious region is introduced outside the two ends as the external boundaries with a width of 2d. The beam is subjected The PD and FE transverse and axial displacements are compared in Fig. 7. There is a very good agreement between PD and FE results. These results verify the accuracy of the current PD formulation for a functionally graded beam subjected to pinned support-roller support boundary condition.

Functionally Graded Beam Subjected to Pinned Support-Pinned Support Boundary Condition
In the second case, the functionally graded beam considered in the previous section is subjected to pinned support at both edges as shown in Fig. 8. Moreover, a horizontal load of P x ¼ 1 N is acting at the center of the beam in addition to the transverse load of P z ¼ 1 N.
PD results for transverse and axial deformations are compared against FE results as shown in Fig. 9 and there is a very good agreement between the two approaches.

Functionally Graded Beam Subjected to Clamped-Clamped Boundary Condition
In the third case, the functionally graded beam is subjected to clamped-clamped boundary condition as shown in Fig. 10. A transverse load of P z ¼ 1 N is applied at the centre of the beam. As demonstrated in Fig. 11, a very good match between PD and FE results is obtained for this particular condition.

Functionally Graded Beam Subjected to Clamped-Free Boundary Condition
In the final numerical case, the functionally graded beam is subjected to clamped-free boundary condition as shown in Fig. 12. A transverse load of P z ¼ 1 N is exerted at the free edge of the beam. As shown in Fig. 13, a very good agreement between PD and FE results is observed which shows that current PD formulation can capture accurate deformation behaviour for functionally graded beams subjected to different types of boundary conditions and loading.

Conclusions
In this study, a new state-based peridynamic formulation was presented for functionally graded Euler beams. The equation of motion was derived by using Lagrange's equations and utilizing Taylor series. To verify the accuracy of the current formulation, four different boundary conditions were considered including pinned support-roller support, pinned support-pinned support, clamped-clamped and clampedfree. For all these boundary conditions, peridynamic results for transverse and axial displacements were compared against finite element analysis results and a very good agreement was obtained between the two approaches. This shows that current peridynamics formulation can capture accurate deformation behaviour for functionally graded beams subjected to different types of boundary conditions and loading. The developed formulation can be further extended to analyze Kirchhoff plates.
Funding Statement: The author(s) received no specific funding for this study.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.