Peridynamic Formulation for Higher-Order Plate Theory

In this study, a new peridynamic model is presented for higher-order plate theory. The formulation is derived by using Euler-Lagrange equation and Taylor’s expansion. The formulation is verified by considering two benchmark problems including simply supported and clamped plates subjected to transverse loading. Moreover, mixed (simply supported-clamped) boundary conditions are also considered to investigate the capability of the current formulation for mixed boundary conditions. Peridynamic results are compared with finite element analysis results and a very good agreement was obtained between the two approaches.


Introduction
Classical continuum mechanics (CCM) developed by Cauchy has been widely used for the analysis of deformation behaviour of materials and structures. Although CCM has been very successful in dealing with numerous complex problems of engineering, it encounters difficulties if the displacement field is discontinuous. This situation mainly arises when cracks occur inside the solution domain. In this case, the partial derivatives in the governing equations of CCM become invalid along the crack surfaces. Moreover, as the technology advances and nanoscale structures become a significant interest, accurate material characteristic at such a small scale cannot always be captured by CCM since CCM does not have a length scale parameter.
To overcome the aforementioned issues related with CCM, a new continuum mechanics formulation, peridynamics, was proposed by Silling [1]. The governing equations of peridynamics (PD) are in the form of integro-differential equations and do not contain any spatial derivatives. Therefore, they are always valid even if the displacement field is discontinuous. Moreover, it has a length scale parameter, horizon, which can be utilized to model structures at nanoscale. According to dell'Isola et al. [2], the origins of peridynamics go back to Piola. Since its introduction, there has been rapid development in peridynamics research especially during the last years. PD has been applied to analyze different material systems including metals [3], composites [4,5], concrete [6] and graphene [7]. Moreover, it is not limited to elasticity behaviour and PD-based plasticity [8], viscoelasticity [9] and viscoplasticity [10] formulations are available. In addition, PD equations have been extended to other fields to perform heat transfer [11], diffusion [12], porous flow [13] and fluid flow [14] analyses. An extensive review on peridynamics is given in Javili et al. [15].
Simplied structures including beams, plates and shells can also be represented in PD framework. Taylor and Steigmann [16] introduced a two-dimensional model for thin plates. Diyaroglu et al. [17] developed a PD Euler-beam formulation which was further extended to Kirchhoff plates by Yang et al. [18]. The effect of transverse shear deformation in thick plates was taken into account by Diyaroglu et al. [19] by developing PD Timoshenko beam and Mindlin plate formulations. O'Grady and Foster [20,21] proposed Euler beam and Kirchhoff plate formulations by utilizing non-ordinary state-based peridynamics. Chowdhury et al. [22] developed a state-based PD formulation suitable for linear elastic shells.
In this study, a new PD formulation is presented for higher-order plate theory suitable for analysis of thick plates. The formulation is developed by using Euler-Lagrange equations and Taylor's expansion. The formulation does not have any limitation on material constants as in bond-based peridynamics. The developed formulation is validated by considering two benchmark problems and peridynamic results are compared against finite element analysis results.

Higher-Order Plate Formulation
The displacement field of any material point, u(x, y, z, t),v(x, y, z, t) and w(x, y, z, t), can be represented in terms of the displacement field of a material point on the mid-plane, u(x, y, 0, t), v(x, y, 0, t) and w(x, y, 0, t), by using Taylor's expansion as where θ x , θ y , θ * x , θ * y and w * are introduced as five new independent variables, respectively, which are defined as (see Fig. 1) In order to simplify the expressions, hereafter, w(x, y, t), θ(x, y, t), w * (x, y, t) and θ * (x, y, t) will be written simply as w, θ, w * and θ * , respectively.
After utilizing the displacement relationships given in Eqs. (2), strain-displacement relationships of 3-dimensional elasticity can be expressed as These strain-displacement relationships can be also expressed in following form by using indicial notation: where i, j = 1, 2, 3 and I, J = 1, 2. Note that this convention where capital letter indices, e.g., I, J, K, … vary from 1 and 2, and lowercase letter indices, e.g., i, j, k, … vary from 1, 2 and 3, will be applied throughout this study. Assuming the material is isotropic and obeys 3-dimensional constitutive relationship, the stress components can be expressed as where C ijkl is the elastic modulus tensor which is defined as with E and ν being elastic (Young's) modulus and Poisson's ratio, respectively. Substituting Eq. (7) into (6) yields: The strain energy density is defined as: Inserting Eqs. (8) and (5) into Eq. (9) and rearranging the indices gives the expression of strain energy density as where ρ is the density, :: u and u represent acceleration and displacement and t denotes time. b(x, t) corresponds to body load vector and is utilized to apply external load. As shown in Fig. 2, f(u ′ , u, x ′ , x) represents the peridynamic interaction (bond) force between two interacting material points located at x and x ′ .
Obtaining analytical solution for Eq. (12) is usually not possible. Instead, numerical techniques such as meshless approach are widely utilized. Therefore, the PD equations of motion for a material point k can be expressed in discrete form as ρ k ð Þ :: where N indicates the total number of family members inside the horizon of material point k and V represents the material point volume.
The PD equations of motion can be derived by utilizing Euler-Lagrange equation as where L = T − U is the Lagrangian. The system's kinetic energy, T, and potential energy, U, can be expressed as and where u (k) and b (k) are the generalized displacement vector and generalized body force density vector, respectively, which in this study can be defined as and Fig. 2 Peridynamic interaction force between two material points [23] b :: Unlike the classical elasticity theory, according to PD theory, the strain energy density function has a non-local characteristic such that the strain energy of a particular material point k depends on both its displacement and all other material points in its family and can be expressed as where u (k) is the displacement vector of 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.
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. (22) and (19) into the Euler-Lagrange equation yields where n 1 = cos φ, n 2 = sin φ with φ being the bond angle with respect to x-axis.
Substituting these PD strain energy expressions into Eq. (23) yields the final PD equations of motion for higher-order plate theory as: Note that when the Poisson's ratio is ν = 0.25, PD equations of motion will have simpler forms:

Boundary Conditions
In this study, the displacement boundary conditions are applied by introducing a fictitious boundary layer, R c , outside the boundary of the actual material domain as shown in Fig. 3. The width of this layer can be chosen as the size of the horizon for the Poisson's ratio value of ν = 1/4 or double size of the horizon for the Poisson's ration of ν ≠ 1/4. This is mainly because the size of the influence domain is equal to double size of the horizon for ν ≠ 1/4. Two common types of boundary conditions, i.e. clamped and simply supported, are explained below for PD higher-order plate formulation.

Clamped Condition
To implement the clamped boundary condition, a fictitious boundary layer is created outside the actual material domain. If the horizon size is chosen as δ = 3Δx, the size of the fictitious domain can be specified as 3Δx for ν = 1/4 or 6Δx for ν ≠ 1/4 in which the discretization size is Δx. The clamped boundary condition constrains zero transverse displacement and zero rotation Fig. 3 Actual, R and fictitious, R c solution domains for the material points 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 θ I and θ * I (I = 1, 2), respectively, to the material points in the fictitious region with respect to the actual displacement field as (see Fig. 4) (i = 1,2,…,6 for ν ≠ 1/4, i = 1,2,3 for ν = 1/4)

Simply Supported Condition
To implement the simply supported boundary condition, a fictitious boundary layer is created outside the actual material domain as in the clamped boundary condition. If the horizon size is chosen as δ = 3Δx, the size of the fictitious domain can be specified as 3Δx for ν = 1/4 or 6Δx for ν ≠ 1/4. The simply supported boundary condition constrains zero transverse displacement for the material points adjacent to the constrained boundary. In this study, this can be achieved by enforcing anti-symmetrical displacement fields for w and w * and symmetrical displacement fields for θ I and θ * I (I = 1, 2), respectively, to the material points in the fictitious region with respect to the actual displacement field as (see Fig. 5)

Numerical Results
In order to verify the PD formulation for a higher-order plate theory, two numerical examples are considered for simply supported and clamped boundary conditions. The PD solutions are compared with the corresponding finite element (FE) analysis results.

Simply Supported Plate Subjected to Transverse Loading
A simply supported plate with a length and width of L = W = 1m and a thickness of h = 0.2m is considered as shown in Fig. 6. The Young's modulus and Poisson's ratio of the plate are E = 200GPa and ν = 1/4, respectively. The model is discretized into one single row of material points along with the thickness direction and the distance between material points is Δx = 1/ 70 m. The horizon size is chosen as δ = 3.015Δx. A fictitious region is introduced outside the edges as the external boundaries with a width of δ. The plate is subjected to a distributed transverse load of p = 100N/m through the y-centre line. The line load is converted to a body load of b ¼ pW 2 w Δx ð ÞΔV ¼ 1:25 Â 10 4 N =m 3 and it is distributed to two columns of material points through the centre line as shown in Fig. 7.
The FE model of the plate is created by using SOLID185 element in ANSYS with 50 elements along the length and width, and 8 elements along the thickness. Boundary conditions below were applied in ANSYS as:  As depicted in Fig. 8, the transverse displacement variation results along the central x-axis and y-axis obtained from PD and FE analyses are compared with each other and a very good agreement is obtained between the two approaches.

Clamped Plate Subjected to Transverse Loading
A clamped plate with a length and width of L = W = 1m and a thickness of h = 0.15m is considered as shown in Fig. 9. Young's modulus and Poisson's ratio of the plate are E = 200GPa and ν = 0.3, respectively. The model is discretized into one single row of material points along with the thickness and the distance between material points is Δx = 1/70 m. The horizon size is chosen as δ = 3.015Δx. A fictitious region is introduced outside the edges as the external boundaries with a width of 2δ. The plate is subjected to a distributed transverse load of p = 100N/m through the y-centre line. The line load is converted to a body load of b ¼ pW 2 w Δx −2 ð Þ ΔV ¼ 1:3021 Â 10 4 N =m 3 and it is distributed to two columns of material points through the centre line as shown in Fig. 10.
The FE model of the plate is created by using the SOLID185 element in ANSYS with 50 elements along the length and width, and 8 elements along the thickness. Boundary conditions below were applied in ANSYS as: PD results for transverse deflection obtained along the central x-axis and y-axis are compared against FE results as shown in Fig. 11. PD results agree very well with FE results.

Plate Subjected to Transverse Loading and Mixed (Simply Supported-Clamped) Boundary Conditions
In the final numerical case, a mixed (simply supported-clamped) boundary condition is considered, as shown in Fig. 12. The plate has a length and width of L = W = 1m and a Transverse deflection results obtained along the central x-axis and y-axis are shown in Fig. 13 and PD and FE results agree very well with each other demonstrating that the current formulation is capable of considering mixed boundary conditions.

Conclusions
In this study, a new peridynamic model was presented for higher-order plate theory. The formulation was derived by using Euler-Lagrange equation and Taylor's expansion. The formulation was verified by considering two benchmark problems including simply supported and clamped plates subjected to transverse loading. Moreover, mixed (simply supportedclamped) boundary conditions are also considered to investigate the capability of the current formulation for mixed boundary conditions. Peridynamic results were compared with finite element analysis results and a very good agreement was obtained between the two approaches. Therefore, it can be concluded that the developed approach can be used as an alternative approach for problems in which higher-order plate theory is applicable.

Appendix
In this section, the derivation of the PD form of strain energy density function given in Eq. (11) is presented. To clarify the derivation, strain energy density expression is separated into four parts: where Transforming W 1 into PD form As shown in Fig. 14 in the Appendix, the function θ can be expressed by using Taylor's expansion up to 1st-order terms about point x as where ξ = |ξ| and unit direction vector n is defined as Multiplying Eq. (A31a) by (A31b) gives Multiplying both sides of Eq. (A33) yields Considering x as a fixed point, integrating both sides of Eq. (A34) over a circular domain with centre of x and radius of δ result in: Transform W 2 into PD form: If we recall the expression given in Eq. (31a) and multiply Eq. (A40) by the direction vector n K yields: Considering x as a fixed point, integrating both sides of Eq. (A41) over a circular domain with centre of x and radius of δ result in: Multiplying both sides of Eq. (A42) with δ IK gives Rewriting Eq. (A43a) with another index, J gives: Following a similar approach, the following expression can also be obtained: 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/.