Analysis of Functionally Graded Timoshenko Beams by Using Peridynamics

In this study, a new peridynamic formulation is presented for functionally graded Timoshenko beams. The governing equations of the peridynamic formulation are obtained by utilising Euler-Lagrange equation and Taylor ’ s expansion. The proposed formulation is validated by considering a Timoshenko beam subjected to different boundary conditions including pinned support-roller support, clamped-roller support and clamped-free boundary conditions. Results from peridynamics are compared against finite element analysis results. A very good agreement is obtained for transverse displacements, rotations and axial displacements along the beam.


Introduction
Peridynamics was introduced by Silling [1] to overcome the limitations of widely used classical continuum mechanics formulation especially for problems including discontinuities in the displacement field due to existence of cracks. Moreover, it has a length scale parameter, horizon, which gives peridynamics a nonlocal characteristic and defines the range of nonlocal interactions between material points. Peridynamics has been utilised to analyse many different types of material systems including metals [2] and composites [3]. Peridynamics can be utilised for not only structural analysis but also for the solution of other physical fields such as heat transfer [4], porous flow [5], diffusion [6], etc. Peridynamics is not limited to elastic behaviour, but peridynamic plasticity [7], viscoelasticity [8] and viscoplasticity [9] formulations are available. Peridynamic formulations are not limited at macro-scale but can be used for the analysis of polycrystalline materials [10] and nano-structures [11]. An extensive review of peridynamics is given in Javili et al. [12].
The original formulation of peridynamics considers only translational degrees of freedom for material points and is capable of performing 3-dimensional analysis. However, this approach can be computationally expensive for certain geometries such as beam, plate and shell-type structures. To capture the correct deformation behaviour of such structures, additional rotational degrees of freedom may be necessary. Such formulations are currently available in the literature. Amongst these Taylor and Steigmann [13] introduced a peridynamic formulation for thin plates. O'Grady and Foster [14,15] developed non-ordinary state-based formulations suitable for Euler-Bernoulli beams and Kirchhoff-Love plate theory. Diyaroglu et al. [16] also proposed a state-based formulation suitable for Euler-Bernoulli beams. Yang et al. [17] extended this formulation to Kirchhoff plates. In order to take into account the transverse shear deformations, Diyaroglu et al. [18] introduced Timoshenko beam and Mindlin plate formulations. In addition, Chowdhury et al. [19] developed a state-based formulation for linear elastic shells.
In this study, a new peridynamic formulation is presented suitable for analysis of functionally graded Timoshenko beams. Functionally graded materials have certain advantages with respect to fibre-reinforced composites since they have continuously varying material properties which does not cause stress concentrations that can lead to delamination failure. There are currently numerous studies available in the literature on analysis of functionally graded Timoshenko beams using classical [20,21] and strain-gradient [22] theories. Peridynamics can be suitable alternative to these approaches. Although there are existing peridynamic formulations available in the literature suitable for the analysis of 2-dimensional and 3-dimensional problems [23,24], the peridynamic functionally graded Timoshenko beam formulation presented in this study has significant benefits in terms of computational time for the analysis of thin and thick functionally graded beams. Proposed peridynamic formulation is derived using Euler-Lagrange equation and Taylor's expansion. To validate the current formulation, several benchmark problems are considered, and peridynamic results are compared against finite element analysis results.

Timoshenko Beam Formulation for Functionally Graded Materials
Timoshenko developed first-order shear deformation theory for thick beams which is also known as Timoshenko beam theory. According to the assumptions of the Timoshenko beam theory, the displacement field of any material point can be expressed in terms of the displacement field of a material point along the central axis in xz plane as (see Fig. 1 where u(x, 0, t) and w(x, 0, t) denote the displacement components of the material point along the central axis in x and z directions, respectively, and t is time. Thus, the strain-displacement relationships can be written as: in which θ ¼ ∂u ∂z represents the rotation angle and u = u(x, 0, t) and w = w(x, 0, t) represent the displacement components of a particular material point located at (x, 0) along the central (neutral) axis of the beam, respectively.
The stress-displacement relationships can also be written as: where E(z) and G(z) denote Young's modulus and shear modulus, respectively, and they are functions of the vertical coordinate, z.
The linear elastic strain energy density of the Timoshenko beam is defined as: where κ is introduced as the shear correction factor. Inserting Eqs. (2) and (3) into (4) yields Peridynamics (PD) is new continuum mechanics formulation introduced by Silling [1]. As opposed to Cauchy's classical continuum mechanics (CCM) formulation, the material points can interact with each other in a nonlocal manner, and the range of interactions is defined by the region called horizon, H. The governing equations of peridynamics are in the form of integro-differential equations and can be written as: where ρ is density, u and :: u represent displacement and acceleration, respectively, and b is the body load vector. t and t ′ define the peridynamic interaction (bond) forces between two material points located at x and x ′ .Analytical solution to Eq. (8) is generally not possible. Numerical techniques are widely used including meshless approach. Therefore, the PD equations of motion for a particular material point k can be expressed in discrete form as: where N indicates the total number of family members and V is the volume of a material point. As shown in Fig. 2, the PD force density vector t (k)(j) represents the force acting on the main material point k by its family member material point j, and, on the contrary, t (j)(k) represents the force acting on the material point j by its family member material point, k.
The PD equations of motion can be derived by utilising Euler-Lagrange equation [25]: where L = T − U is the Lagrangian which is the difference between the system's kinetic energy, T, and potential energy, U, that can be expressed as: and with u and b being the generalised displacement vector and generalised body force density vector, respectively. In this study, they can be defined as: Here, the entries of the body force density vector, b x , b θ and b z , correspond to axial force along x-axis, bending moment and transverse force densities, respectively. Considering rectangular cross-section case, inserting Eq. (1) into (11), integrating throughout the thickness and dividing by the thickness give the kinetic energy of the system as: The first term of Euler-Lagrange equation can be obtained by substituting Eq. (14) into Eq. (10) as: Unlike the classical elasticity theory, according to PD theory, the strain energy density function has nonlocal form such that the strain energy of a particular material point k depends on both its displacement and displacements of all other material points in its horizon, which can be expressed as: where u (k) is the displacement vector of the material point k and u i k ð Þ (i = 1, 2, 3, ⋯) is the displacement vector of the ith material point inside the horizon of the material point k. Therefore, 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: Thus, the second term of the Euler-Lagrange equation can be written as: Inserting Eqs. (15) and (18) The strain energy density function given in Eq. (7) can be written in nonlocal PD form for a particular material point k. To achieve this, it is necessary to transform all differential terms into equivalent PD forms in accordance with the form of PD strain energy density given in Eq. (17). As derived in Appendix 1, the strain energy density function of the material point k and its family member j can be expressed as: Substituting Eqs. (20a) and (20b) into (18) and renaming the summation indices, the second term of Euler-Lagrange equation becomes: Therefore, the entire set of equations of motion is given by inserting Eqs. (15) and (21) ::

Application of Boundary Conditions
Application of boundary conditions is a critical step of peridynamic analysis. Since the PD equations of motion are in the form of integro-differential equations, the approach to apply boundary conditions is different than the approach used in CCM. In this study, to apply displacement boundary conditions, a fictitious region, R c , is created outside of the actual solution domain, R. The width of the fictitious layer can be chosen as the same size of the horizon. Two common types of boundary conditions including clamped and simply supported boundary conditions are presented below for PD Timoshenko beam formulation for functionally graded materials.

Clamped Boundary Condition
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Δx in which the discretisation size is Δx. Hence, the size of the fictitious region is chosen as δ.
In Timoshenko beam theory, the clamped boundary condition imposes zero transverse displacement, zero axial displacement and zero rotation for the material point adjacent to the clamped end as: In our study, these conditions can be achieved by enforcing mirror image of the transverse displacement field for the material points adjacent to the clamped end and anti-symmetric image of rotational and axial displacements fields, as shown in Fig. 3, such that: Note that enforcing conditions given in Eqs. (24a, b) also implicitly satisfies the conditions related with accelerations of the material point adjacent to the clamped end, i.e.

Simply Supported Boundary Condition
To implement the simply supported boundary condition, the fictitious layer is introduced outside the real material domain, whose size is again chosen to be equal to δ. In Timoshenko beam theory, the simply supported boundary condition imposes zero transverse displacement and zero curvature for the material point adjacent to the clamped end as: In our study, the above conditions can be satisfied by enforcing anti-symmetrical transverse displacement field and symmetrical rotational displacement field to the material points in the fictitious region with respect to the actual displacement field, as shown in Fig. 4, which is defined as: Fig. 9 Material variation in the thickness direction for the finite element model The discussion about simply supported boundary condition for axial deformation should be separated into two parts, which are pinned support and roller support, respectively. The implementation of axial deformation for pinned support is achieved by enforcing antisymmetrical displacement field for fictitious material domain with respect to actual material domain, as shown in Fig. 4, which can be defined as: On the other hand, the implementation for roller support boundary condition requires symmetrical displacement field adjacent to the boundary as (see Fig. 5)

Numerical Results
To verify the validity of the PD formulation for functionally graded Timoshenko beams, the PD solutions are compared with the corresponding finite element (FE) analysis results. In this study, the functionally graded material properties are chosen as Young's modulus E(z) and shear modulus G(z), and they are 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 denotes the total thickness of the beam as shown in Fig. 6. All numerical examples considered

Timoshenko Beam with Pinned Support-Roller Support Boundary Conditions
A simply supported functionally graded beam with length, width and thickness of L = 1 m, W = 0.01 m and t = 0.1 m is considered as shown in Fig. 7. The beam is constrained by pinned The beam is subjected to a concentrated transverse load of P z = 100 N at the central point. The load is converted to a body load of b ¼ Pz 2ΔV ¼ 2:5e7 N=m 3 , and it is applied to two central material volumes of the model (see Fig. 8).
The FE model of the beam is created by using SHELL181 element in ANSYS with dimensions of 1 m × 0.01 m × 0.1 m. The FE model is meshed with 100 elements along the length. To model the functionally graded beam, the model is divided to 50 layers with varying homogeneous materials properties through the thickness. Young's modulus varies linearly The PD and FE transverse displacements, w(x), rotations, θ(x), and axial displacements, u(x), are compared in Fig. 10. Both approaches yield similar displacement and rotation variations. These results verify the accurateness of the current PD formulation for a functionally graded beam theory for pinned support-roller support boundary condition.

Timoshenko Beam with Clamped-Roller Support Boundary Conditions
The Timoshenko beam with clamped-roller supported boundary conditions is under consideration (see Fig. 11). The beam has the same geometrical and elastic properties as previous case. The model is discretised into one single row of material points along with the thickness, and the distance between material points is Δx = 0.002 m. A fictitious region is introduced outside the ends as the external boundaries with a width of δ. The beam is subjected to a concentrated transverse load of P z = 100 N at the central point. The load is converted to a body load of b ¼ P z 2ΔV ¼ 2:5e7 N =m 3 , and it is subjected to two central material volumes of the model (see Fig. 12).
As depicted in Fig. 13, for transverse displacements, w(x), rotations, θ(x), and axial displacements, u(x), a very good agreement is obtained between PD and FE results.

Timoshenko Beam with Clamped-Free Boundary Conditions
The performance of a cantilever beam is under investigation (see Fig. 14). The beam has the same geometrical and elastic properties as previous case. The model is discretised into one single row of material points along with the thickness, and the distance between material points is Δx = 0.002 m. Fictitious regions are introduced outside the ends as the external boundaries with a width of δ. The fictitious region at the free end is included to improve the accuracy of the calculation of the material points adjacent to the free end due to lack of interactions at the free end. The beam is subjected to a bending moment of M = 100 Nm at the free end. The load is converted to a body load of b θ ¼ M ΔV ¼ 5e7 N=m 2 , and it is applied to the last material point adjacent to the free end (see Fig. 15).
The transverse displacements, w(x), rotations, θ(x), and axial displacements, u(x), along the beam are obtained from both PD and FE analysis. As shown in Fig. 16, a very good agreement is observed between the two approaches.

Conclusions
In this study, a new peridynamic formulation was presented for functionally graded Timoshenko beams. The proposed formulation was validated by considering a Timoshenko beam subjected to different boundary conditions including pinned support-roller support, clampedroller support and clamped-free boundary conditions. Results from peridynamics were compared against finite element analysis results. A very good agreement was obtained for transverse displacements, rotations and axial displacements along the beam. Therefore, it can be concluded that peridynamics can be a suitable alternative for the analysis of functionally graded Timoshenko beams.
In order to obtain the strain energy density function in PD form, it is necessary to transform each term into an equivalent PD expression, and this can be achieved by utilising Taylor's expansion. If we expand the axial displacement u about a point x and ignore the higher order terms: By squaring both sides of Eq. (31) and then dividing each term by |ξ| yields: Considering x as a fix point and integrating each term of Eq. (32) over the region (−δ, δ) results in: which then gives: Equation (33b) can be written in discretised form as: where A represents the cross-sectional area of the beam. Two other terms can be obtained similarly as: Similarly, PD form of strain energy density function for material point j is: