Peridynamic Mindlin Plate Formulation for Functionally Graded Materials

: In this study, a new peridynamic Mindlin plate formulation is presented which is suitable for the analysis of functionally graded materials. The governing equations of peridynamic formulation are obtained by using Euler ‐ Lagrange equations in conjunction with Taylor’s expansion. To validate the new formulation, three different numerical benchmark problems are considered for a Mindlin plate subjected to simply supported, fully clamped and mixed (clamped ‐ simply supported) boundary conditions. Peridynamic results are compared against results from finite element analysis and a good agreement is observed between the two methods.


Introduction
As the manufacturing technology advances, it becomes possible to design materials that have better properties with respect to traditional materials including metals and fibre-reinforced composite materials. Usage of fibre-reinforced composite materials is increasing in different industries including automotive, aerospace and marine fields. Although fibre-reinforced composite materials are light, have good impact properties and corrosion resistance, they suffer from delamination damage occurring between neighbouring plies which significantly reduces the load carrying capacities of these materials. Delamination damage can be avoided by continuously varying material properties in the form of functionally graded materials.
There are currently numerous studies in the literature focusing on plate formulations of functionally graded materials. Amongst these Vel and Batra [1] provided three-dimensional exact solution for the vibration of functionally graded rectangular plates. Shen [2] obtained nonlinear bending response of functionally graded plates subjected to transverse loading by using Reddy's higher order shear deformation plate theory. Zenkour [3] presented generalised shear deformation theory for bending analysis of functionally graded plates. Bian et. al. [4] derived analytical solutions for functionally graded plates under cylindrical bending using first-and third-order shear deformation theories. Carrera et. al. [5] investigated the effect of thickness stretching in functionally graded plates by using Carrera's Unified Formulation. Ferreira et. al. [6] used meshless method and third-order shear deformation theory to analyse functionally graded plates. Kashtalyan [7] utilised Plevako general solution of the equilibrium equations for inhomogeneous isotropic media to obtain three-dimensional elasticity solution for functionally graded simply supported plates. Xiang and Kang [8] performed bending analysis of functionally graded plates by using nth-order shear deformation theory and meshless global collocation method based on the thin plate spline radial basis function. Thus, the strain-displacement relationship can be written as which can be also expressed in terms of indicial notation as where κ is introduced as shear coefficient. Note that the subscript indices, I, J, ⋯ 1 ,2 , and this convention will be applied throughout this study.
For planar isotropic materials, the stress-strain relationships can be written as: where is the shear modulus.
Note that the transverse normal stress, , is considered to be small compared to in-plane stresses. Thus, it is discarded from the stress components set and this simplifies the 3-dimensional Hooke's law and makes it a 2-dimensional plane-stress material constitutive law. The stress components can also be expressed in indicial notation as:

Peridynamic Mindlin Plate Formulation
Peridynamics (PD) is a non-local continuum mechanics formulation. Therefore, material points inside the solution domain can interact with each other material points in a non-local manner. The range of non-local interactions is defined as 'horizon', . The equations of motion of peridynamics can be written as or in discrete form for the material point as (10) where and are the peridynamic force densities between two material points k and j as shown in Figure 2 and N is the number of material points inside the horizon. Note that the PD equations of motion given in Equations (9) and (10) have the same form as in classical continuum mechanics except the integration or summation term.
where is the Lagrangian. The kinetic energy per unit area, , can be expressed as where is the generalised displacement vector, which in this study can be defined as (13) Substituting Equation (1) into Equation (12) results in The total kinetic energy of the body can be casted by integrating Equation (14) over the whole mid-plane as which can be written in discretised form as where is the area that contains the material point . The PD strain energy density function has a non-local form such that the strain energy of a certain material point depends on both its displacement and all other material points in its family, which can be expressed as where is the displacement vector of material point and ( 1, 2, 3, ⋯ ) is the displacement vector of the th material point within the horizon of the material point . Similar to Equation (15), 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 where is the body force density vector, which in this study has the following components (18) with , and correspond to in-plane loads, moments and transverse body loads, respectively. The first term of the Euler-Lagrange's equation can be obtained by inserting Equations (15b) and (17) into Equation (11) as Similarly, the second term of the Euler-Lagrange's equation becomes Inserting Equations (19) and (20) into the Euler-Lagrange's equation yields: In order to write the non-local form of strain energy function of the material point , Equation (8b), it is necessary to transform all the local terms into an equivalent PD form by also considering PD strain energy expression given in Equation (16). As derived in Appendix A, the strain energy density function of the material point and its family member can be expressed as cos , sin and is the orientation of the peridynamic interaction with respect to horizontal axis.
Inserting Equations (22a) and (22b) into Equation (21) results in complete PD equations of motion for functionally graded Mindlin plates as In particular, when Poisson's ratio, ν z , Equation (23) can be reduced to PD bond-based formulation as

Numerical Results
To verify the validity of the PD formulation for functionally graded Mindlin plates, 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, and shear modulus and they are assumed to vary linearly through the thickness as where and denote the Young's modulus of the top and bottom surfaces of the plate, and ℎ denotes the total thickness of the plate. The shear correction coefficient is chosen as .
In the following three numerical cases, a square plate with length and width of 1m and thickness of ℎ 0.15m is considered. The plate is subjected to different boundary conditions. The Young's modulus of the top and bottom surface are chosen as 200 and 100 . The PD models are discretized into one single row of material points through the thickness and 65 65 material points throughout the xy plane. Thus, the distance between two adjacent material points is ∆x and the area attached on each material points is ∆A ∆ . A fictitious region is introduced outside the edges as the external boundaries with a width of 6∆x to apply boundary conditions as explained in Appendix B. The horizon size can be approximately chosen as 3∆ .
The corresponding FE models are created in ANSYS by using SHELL181 elements with 50 50 elements throughout the plate. In order to obtain the functionally graded character, the model is divided into 50 layers with varying homogeneous material properties throughout the thickness. The Young's modulus varies gradually over the thickness from the first layer 101 to the last layer 199 , as shown in Figure 3. The Poisson's ratio, ν 0.3, is applied in ANSYS.

Simply Supported Functionally Graded Mindlin Plate
In the first example case, a simply supported functionally graded Mindlin plate subjected to a distributed load of 100,000N/m through the central line is taken into consideration (see Figure   4). For the PD model, the load is transformed into a body load of   Assuming the coordinate system is attached at the geometric centre of the plate, the following boundary conditions are applied in ANSYS: The PD results for in-plane and transverse displacements, and rotations are obtained and compared with FEA results for the material points located along central x-and y-axes. As depicted in Figure 6, PD and FEA results agree very well with each other.

Fully Clamped Functionally Graded Mindlin Plate
In the second example case, the functionally graded Mindlin plate considered in the previous example is subjected to fully clamped boundary condition (see Figure 7). In ANSYS, the clamped boundary condition is achieved by constraining all degrees of freedom along the external boundaries. Based on the comparison between the PD and FEA results as shown in Figure 8, it can be concluded that current PD formulation can also provide accurate results for fully clamped boundary condition.

Functionally Graded Mindlin Plate Subjected to Mixed Boundary Conditions
The last numerical case aims to verify the current PD formulation for mixed boundary conditions, i.e., clamped-simply supported. As shown in Figure 9, edges along the horizontal direction are subjected to clamped boundary conditions whereas the remaining two edges are subjected to simply supported boundary conditions. Assuming the coordinate system is attached at the geometric centre of the plate, the following boundary conditions are applied in ANSYS: As depicted in Figure 10, also for this mixed-boundary conditions case, PD and FEM results agree well with each other.

Conclusions
In this study, a new peridynamic formulation was presented for functionally graded Mindlin plates. The governing equations were obtained by using Euler-Lagrange equations. To validate the formulation, three different benchmark problems were considered for simply supported, fully clamped and mixed (clamped-simply supported) boundary conditions. Peridynamic results were compared against finite element analysis results and a good agreement is observed for both in-plane and transverse displacements, and rotations. Hence, it can be concluded that the present peridynamic formulation can be a suitable alternative for the analysis of functionally graded Mindlin plates subjected to different types of boundary conditions. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
As explained above, according to classical continuum mechanics, the strain energy per unit area of functionally graded Mindlin plates can be written as In order to obtain the strain energy function in PD form, it is necessary to transform each local term in Equation (A1) into their equivalent non-local form. This can be achieved by using Taylor expansion.
and added with (A17) yields: Rewriting Equation (A20a) with a different index results in: Multiplying Equation (A20a) with (A20b) and then dividing each term by ξ yields:  To implement the clamped boundary condition, a fictitious boundary layer is created outside the actual material domain (see Figure A3). The horizon size can be approximately chosen as 3∆ in which the discretisation size is ∆ . The size of the fictitious region can be specified to be equal to 6∆x for 1/3 and 3∆x for 1/3.
According to classical theory, the clamped boundary condition imposes zero in-plane displacements, zero transverse displacement and zero rotations for the material point adjacent to the clamped end as In this 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 in-plane and rotational displacements fields as following * To implement the simply supported boundary condition, the fictitious layer is introduced outside the actual material domain (see Figure A4), whose size is again chosen to be equal to 6∆x for 1/3 and 3∆x for 1/3. According to classical theory, the simply supported boundary condition can be achieved by imposing zero in-plane displacements, zero transverse displacements and zero curvature for the material point adjacent to the boundaries as