A state-based peridynamic formulation for functionally graded Kirchhoff plates

Abstract

Functionally graded materials are a potential alternative to traditional fibre-reinforced composite materials as they have continuously varying material properties which do not cause stress concentrations. In this study, a state-based peridynamic model is presented for functionally graded Kirchhoff plates. Equations of motion of the new formulation are obtained using the Euler–Lagrange equation and Taylor’s expansion. The formulation is verified by considering several benchmark problems including a clamped plate subjected to transverse loading and a simply supported plate subjected to transverse loading and inclined loading. The material properties are chosen such that Young’s modulus is assumed to be varied linearly through the thickness direction and Poisson’s ratio is constant. Peridynamic results are compared against finite element analysis results, and a very good agreement is obtained between the two approaches.

1. Introduction

Functionally graded materials (FGMs) are a potential alternative to traditional fibre-reinforced composite materials (FRCMs) as they have continuously varying material properties which do not cause stress concentrations. Such stress concentrations arise in FRCMs and cause delamination-type failure. With the advancement of additive manufacturing technologies, it is expected that the usage of FGMs will increase in the near future.

There are many studies in the literature focusing on analysis of FGMs. Amongst these Vel and Batra [1] presented three-dimensional exact solutions for free and forced vibrations of simply supported functionally graded plates using a power series method and considering classical plate theory, first-order shear deformation theory and third-order shear deformation theory. Zenkour [2] studied the bending deformation of a simply supported functionally graded plate subjected to transverse uniform load by using generalised shear deformation theory without using shear correction factors. Shen [3] performed nonlinear bending analysis for simply supported functionally graded plates subjected to transverse loading by taking into account the effect of temperature and using Reddy’s higher-order shear deformation plate theory. Bian et al. [4] utilised a plate theory using the concept of shape function of the transverse coordinate parameter to investigate functionally graded plates subjected to cylindrical bending. Carrera et al. [5] examined the effect of thickness stretching in functionally graded plates by utilising Carrera’s Unified Formulation. As an alternative formulation, a new continuum mechanics formulation, peridynamics (PD), can be utilised.

PD was introduced by Silling [6] to overcome the limitations of Classical Continuum Mechanics. Governing equations of PD is in the form of integro-differential equations and does not contain spatial derivatives. Therefore, these equations are still valid even if the displacement field is discontinuous due to existence of cracks. Moreover, PD is a non-local continuum mechanics formulation and it has a length-scale parameter, horizon, which defines range of non-local interactions between material points. PD has been used for analysing different types of materials including metals and composites [7–10] and is suitable for multiphysics analysis [11–15]. Moreover, several peridynamic formulations are available in the literature to analyse beams, plates and shells [16–21].

In this study, PD is utilised to investigate the bending behaviour of functionally graded plates based on Kirchhoff plate theory for the first time in the literature. The governing equations are obtained by using Euler–Lagrange equation and Taylor’s expansion. Several benchmark problems are considered to validate the newly developed formulation.

2. Functionally graded Kirchhoff plate

Kirchhoff developed a complete and adequate set of equation for linear theory of thin plates. According to the Kirchhoff plate theory, it is assumed that a transverse normal to the mid-surface in the undeformed state remains straight, normal to the mid-surface and it has no change in length during deformation. Moreover, the displacement field of any material point can be represented in terms of the displacement field of the material points at the mid-surface as

$$u_1(x_1,x_2,z)={\overline{u}}_1(x_1,x_2)+z·{\overline{u}}_1(x_1,x_2{)}_{,z}$$

(1a)

$$u_2(x_1,x_2,z)={\overline{u}}_2(x_1,x_2)+z·{\overline{u}}_2(x_1,x_2{)}_{,z}$$

(1b)

$$w(x_1,x_2,z)=\overline{w}(x_1,x_2)$$

(1c)

where ${\overline{u}}_1$ , ${\overline{u}}_2$ and $\overline{w}$ denote the displacement of a material point at the mid-surface in $x_1($ x), $x_2(y$ ) and z directions, respectively. Thus, the strain-displacement relationships can be written as

$${\epsilon }_{11}=\frac{\partial {\overline{u}}_1}{\partial x_1}+z\frac{\partial {\theta }_1}{\partial x_1}$$

(2a)

$${\epsilon }_{22}=\frac{\partial {\overline{u}}_2}{\partial x_2}+z\frac{\partial {\theta }_2}{\partial x_2}$$

(2b)

$${\epsilon }_{33}=0$$

(2c)

$${\epsilon }_{12}={\epsilon }_{21}=\frac{1}{2}[\frac{\partial {\overline{u}}_1}{\partial x_2}+\frac{\partial {\overline{u}}_2}{\partial x_1}+z(\frac{\partial {\theta }_1}{\partial x_2}+\frac{\partial {\theta }_2}{\partial x_1})]$$

(2d)

$${\epsilon }_{13}={\epsilon }_{31}=\frac{1}{2}({\theta }_1+\frac{\partial \overline{w}}{\partial x_1})=0$$

(2e)

$${\epsilon }_{23}={\epsilon }_{32}=\frac{1}{2}({\theta }_2+\frac{\partial \overline{w}}{\partial x_2})=0$$

(2f)

in which ${\theta }_1=\frac{\partial {\overline{u}}_1}{\partial z}$ and ${\theta }_2=\frac{\partial {\overline{u}}_2}{\partial z}$ represent the rotation angles of a particular material point at the mid-surface which can be expressed in terms of $\overline{w}$ from equations (2e) and (2f) as

$${\theta }_1=-\frac{\partial \overline{w}}{\partial x_1}$$

(3a)

and

$${\theta }_2=-\frac{\partial \overline{w}}{\partial x_2}$$

(3b)

Therefore, eliminating ${\theta }_1$ and ${\theta }_2$ , strain components in a thin plate can be easily represented in terms of only three mid-surface parameters, ${\overline{u}}_1$ , ${\overline{u}}_2$ and $\overline{w}$ , as

$$\mathbf{\epsilon }=\left[\begin{array}{cc}{\epsilon }_{11}\hfill & {\epsilon }_{12}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{21}}& {\epsilon }_{22}\hfill \end{array}\right]=\left[\begin{array}{cc}{\overline{u}}_{1,1}-z·{\overline{w}}_{,11}\hfill & \frac{1}{2}({\overline{u}}_{1,2}+{\overline{u}}_{2,1})-z·{\overline{w}}_{,12}\hfill \\ \multicolumn{1}{c}{\frac{1}{2}({\overline{u}}_{1,2}+{\overline{u}}_{2,1})-z·{\overline{w}}_{,12}}& {\overline{u}}_{2,2}-z·{\overline{w}}_{,22}\hfill \end{array}\right]$$

(4)

or

$${\epsilon }_{\mathrm{ij}}=\frac{1}{2}({\overline{u}}_{i,j}+{\overline{u}}_{j,i})-z·{\overline{w}}_{,\mathrm{ij}}$$

(5)

where the indices $i$ and $j$ vary from 1 to 2. Note that from this stage, for the sake of conciseness, the notation ‘’overbar”, $\overline{}$ , will be removed from displacement components.

Plane stress constitutive relation is held in the thin plate of functionally graded material as

$$\mathbf{\sigma }=\mathit{C}\mathbf{\epsilon }$$

(6a)

and

$$\left\{\begin{array}{c}{\sigma }_{11}\hfill \\ \multicolumn{1}{c}{{\sigma }_{22}}\\ \multicolumn{1}{c}{{\sigma }_{12}}\end{array}\right\}=\frac{E(z)}{1-\nu {(z)}^2}\left[\begin{array}{ccc}1\hfill & \nu (z)\hfill & 0\hfill \\ \multicolumn{1}{c}{\nu (z)}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1-\nu (z)}{2}\hfill \end{array}\right]\left\{\begin{array}{c}{\epsilon }_{11}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{22}}\\ \multicolumn{1}{c}{2{\epsilon }_{12}}\end{array}\right\}$$

(6b)

where $E(z)$ and $\nu (z)$ represent elastic modulus and Poisson’s ratio, respectively, and both of them vary in the thickness direction z. Equation (6b) can also be written in indicial notation as

$${\sigma }_{\mathrm{ij}}=C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}$$

(6c)

where $C_{\mathrm{ijkl}}$ is the components of the material property tensor and it is defined as

$$C_{\mathrm{ijkl}}=\frac{E(z)}{1-\nu {\left(z\right)}^2}[\frac{1-\nu \left(z\right)}{2}({\delta }_{\mathrm{il}}{\delta }_{\mathrm{jk}}+{\delta }_{\mathrm{ik}}{\delta }_{\mathrm{jl}})+\nu \left(z\right){\delta }_{\mathrm{ij}}{\delta }_{\mathrm{kl}}]$$

(7)

The linear elastic strain energy density of the plate can be expressed as

$$W=\frac{1}{2}{\sigma }_{\mathrm{ij}}{\epsilon }_{\mathrm{ji}}=\frac{1}{2}{C}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}{\epsilon }_{\mathrm{ji}}$$

(8)

Substituting equations (6) and (7) into equation (8) and rearranging the terms results in

$$W=\frac{E(z)}{2(1-\nu {\left(z\right)}^2)}\left\{\begin{array}{c}\frac{1-\nu \left(z\right)}{2}(u_{i,j}{u}_{i,j}+u_{i,j}{u}_{j,i}+u_{i,i}{u}_{j,j})+\frac{3\nu \left(z\right)-1}{2}{u}_{i,i}{u}_{j,j}\hfill \\ \multicolumn{1}{c}{+z^2(w_{,\mathrm{ii}}{w}_{,\mathrm{jj}}+(1-\nu \left(z\right))(w_{,\mathrm{ij}}{w}_{,\mathrm{ij}}-w_{,\mathrm{ii}}{w}_{,\mathrm{jj}}))}\\ \multicolumn{1}{c}{-2z[u_{i,i}{w}_{,\mathrm{jj}}+(1-\nu \left(z\right))(\frac{1}{2}(u_{i,j}+u_{j,i})w_{,\mathrm{ij}}-u_{i,i}{w}_{,\mathrm{jj}})]}\end{array}\right\}$$

(9)

As the strain energy density changes in the thickness direction, the average strain energy density can be obtained by integrating equation (9) over the thickness and dividing by thickness, which yields

$$W=\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{2(1-\nu {\left(z\right)}^2)}\left\{\begin{array}{c}\frac{1-\nu \left(z\right)}{2}(u_{i,j}{u}_{i,j}+u_{i,j}{u}_{j,i}+u_{i,i}{u}_{j,j})+\frac{3\nu \left(z\right)-1}{2}{u}_{i,i}{u}_{j,j}\hfill \\ \multicolumn{1}{c}{+z^2(w_{,\mathrm{ii}}{w}_{,\mathrm{jj}}+(1-\nu \left(z\right))(w_{,\mathrm{ij}}{w}_{,\mathrm{ij}}-w_{,\mathrm{ii}}{w}_{,\mathrm{jj}}))}\\ \multicolumn{1}{c}{-2z[u_{i,i}{w}_{,\mathrm{jj}}+(1-\nu \left(z\right))(\frac{1}{2}(u_{i,j}+u_{j,i})w_{,\mathrm{ij}}-u_{i,i}{w}_{,\mathrm{jj}})]}\end{array}\right\}\mathrm{dz}$$

(10)

Note that equation (10) is composed of three independent parts. The first and second rows represent the strain energy densities that occur due to in-plane and flexural (bending) deformation, respectively, and the last row arises due to coupling of in-plane and flexural deformations. If the material properties are constant over the thickness, the last row of equation (10) will be cancelled out after integration and equation (10) will become uncoupled. Thus, the strain energy density of the plate can be written concisely as

$$W=W_I+W_{\mathrm{II}}+W_{\mathrm{III}}$$

(11)

where $W_I$ , $W_{\mathrm{II}}$ and $W_{\mathrm{III}}$ denote the strain energy densities due to in-plane deformation, flexural deformation and coupled deformation, respectively, which can be written as

$$W_I=\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{2(1-\nu {\left(z\right)}^2)}[\frac{1-\nu \left(z\right)}{2}(u_{i,j}{u}_{i,j}+u_{i,j}{u}_{j,i}+u_{i,i}{u}_{j,j})+\frac{3\nu \left(z\right)-1}{2}{u}_{i,i}{u}_{j,j}]\mathrm{dz}$$

(12a)

$$W_{\mathrm{II}}=\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z^2}{2(1-\nu {\left(z\right)}^2)}[w_{,\mathrm{ii}}{w}_{,\mathrm{jj}}+(1-\nu \left(z\right))(w_{,\mathrm{ij}}{w}_{,\mathrm{ij}}-w_{,\mathrm{ii}}{w}_{,\mathrm{jj}})]\mathrm{dz}$$

(12b)

$$W_{\mathrm{III}}=-\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z}{1-\nu {(z)}^2}\{u_{i,i}{w}_{,\mathrm{jj}}+(1-\nu \left(z\right))(\frac{1}{2}(u_{i,j}+u_{j,i})w_{,\mathrm{ij}}-u_{i,i}{w}_{,\mathrm{jj}})\}\mathrm{dz}$$

(12c)

3. Peridynamic formulation for functionally graded Kirchhoff plate

PD is a non-local theory, which is different than the classical local continuum theory, as the state of each material point is not only influenced by the material points located in its immediate vicinity but is also influenced by material points which are located within a region of a finite radius named as ‘horizon’, $H$ . The equation of motion in PD can be written as

$$\rho \left(\mathit{x}\right)\stackrel{··}{\mathit{u}}(\mathit{x},t)={\int }_H\mathbf{f}(\mathit{u}\prime -\mathit{u},\mathit{x}\prime -\mathit{x})dV\prime +\mathit{b}(\mathit{x},t)$$

(13)

where $\rho$ and t represent density and time, respectively, $\mathit{u}$ and $\stackrel{··}{\mathit{u}}$ denote displacement and $\mathit{b}$ is the body load vector. In equation (13), $\mathbf{f}$ represents the interaction (bond) force vector between material points located at $\mathit{x}$ and $\mathit{x}\prime$ . The closed-form solution to equation (13) is generally not available. Numerical approaches such as the meshless approach are widely used to solve equation (13). Therefore, the PD equations of motion for a particular material point $k$ can be expressed in summation form as

$${\rho }_{(k)}{\stackrel{··}{\mathit{u}}}_{(k)}={\sum }_{j=1}^N{\mathit{f}}_{(k)(j)}{V}_{(j)}+{\mathit{b}}_{(k)}$$

(14)

where $N$ indicates the total number of material points inside the horizon of the material point, k. The interaction force vector, ${\mathit{f}}_{(k)(j)}$ , between material points $k$ and $j$ is defined as

$${\mathit{f}}_{(k)(j)}={\mathit{t}}_{(k)(j)}-{\mathit{t}}_{(j)(k)}$$

(15)

As shown in Figure 1, the PD force density vector ${\mathit{t}}_{(k)(j)}$ represents the force acting on the main material point $k$ by its family member material point $j$ , and, on the contrary, ${\mathit{t}}_{(j)(k)}$ represents the force acting on material point $j$ by its family member material point, $k$ .

Figure 1.

Peridynamic interaction force between material points [22].

Unlike the classical elasticity theory, according to PD theory, the strain energy density function of a particular material point $k$ depends on the relative displacement between $k$ and all other material points in its horizon which can be expressed as

$$W_{(k)}=W_{\left(k\right)}({\mathit{u}}_{\left(1^k\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{u}}_{\left(2^k\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{u}}_{\left(3^k\right)}-{\mathit{u}}_{\left(k\right)},\cdots )$$

(16)

where ${\mathit{u}}_{\left(k\right)}$ is the displacement vector of point $k$ and ${\mathit{u}}_{\left(i^k\right)}$ ( $i=1,2,3,\cdots$ ) is the displacement vector of the $i$ th material point within 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

$$U={\sum }_k{W}_{\left(k\right)}({\mathit{u}}_{\left(1^k\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{u}}_{\left(2^k\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{u}}_{\left(3^k\right)}-{\mathit{u}}_{\left(k\right)},\cdots )V_{(k)}-\sum _k{\mathit{b}}_{(k)}{\mathit{u}}_{(k)}{V}_{(k)}$$

(17)

where ${\mathit{b}}_{(k)}$ is the body load vector of the material point $k$ with a unit of ‘force per unit volume’. Similarly, the total kinetic energy of the body can be obtained by summing kinetic energies of all material points as

$$T=\frac{1}{2}{\sum }_k{\rho }_{(k)}{\stackrel{·}{\mathit{u}}}_{(k)}·{\stackrel{·}{\mathit{u}}}_{(k)}{V}_{(k)}$$

(18)

Thus, the Lagrangian of the body can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{L=T-U=\frac{1}{2}{\sum }_k{\rho }_{\left(k\right)}{\stackrel{·}{\mathit{u}}}_{\left(k\right)}·{\stackrel{·}{\mathit{u}}}_{\left(k\right)}{V}_{\left(k\right)}-{\sum }_k{W}_{\left(k\right)}({\mathit{u}}_{\left(1^k\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{u}}_{\left(2^k\right)}-{\mathit{u}}_{\left(k\right)},{\mathit{u}}_{\left(3^k\right)}-{\mathit{u}}_{\left(k\right)},\cdots )V_{\left(k\right)}}\\ +{\sum }_k{\mathit{b}}_{(k)}{\mathit{u}}_{(k)}{V}_{(k)}\hfill \end{array}$$

(19)

Utilising the Lagrangian term given in equation (19), the equations of motion can be obtained from Euler–Lagrange equation as

$$\frac{d}{\mathrm{dt}}\frac{\partial L}{\partial {\stackrel{·}{\mathit{u}}}_{\left(k\right)}}-\frac{\partial L}{\partial {\mathit{u}}_{\left(k\right)}}=0$$

(20)

which yields

$${\rho }_{\left(k\right)}{\stackrel{··}{\mathit{u}}}_{\left(k\right)}-[{\sum }_{i^k}\frac{\partial W_{\left(k\right)}}{\partial ({\mathit{u}}_{\left(i^k\right)}-{\mathit{u}}_{\left(k\right)})}\frac{\partial ({\mathit{u}}_{\left(i^k\right)}-{\mathit{u}}_{\left(k\right)})}{\partial {\mathit{u}}_{\left(k\right)}}+\sum _{i^k}{\sum }_{i^{k^k}}\frac{\partial W_{\left(i^k\right)}}{\partial ({\mathit{u}}_{\left(i^{k^k}\right)}-{\mathit{u}}_{\left(i^k\right)})}\frac{\partial ({\mathit{u}}_{\left(i^{k^k}\right)}-{\mathit{u}}_{\left(i^k\right)})}{\partial {\mathit{u}}_{\left(k\right)}}\frac{V_{\left(i^k\right)}}{V_{\left(k\right)}}+{\mathit{b}}_{(k)}]=0$$

(21)

where point $i^k$ is a family member of the material point $k$ and similarly point $i^{k^k}$ is a family member of the material point $i^k$ . Note that the following relationship holds

$$\frac{\partial ({\mathit{u}}_{\left(i^{k^k}\right)}-{\mathit{u}}_{\left(i^k\right)})}{\partial {\mathit{u}}_{\left(k\right)}}={\delta }_{i^{k^k}k}$$

(22)

Substituting equation (22) into equation (21) and rearranging the summation index leads to

$${\rho }_{\left(k\right)}{\stackrel{··}{\mathit{u}}}_{\left(k\right)}={\sum }_j\frac{\partial W_{\left(k\right)}}{\partial ({\mathit{u}}_{\left(j\right)}-{\mathit{u}}_{\left(k\right)})}-\sum _j\frac{\partial W_{\left(j\right)}}{\partial ({\mathit{u}}_{\left(k\right)}-{\mathit{u}}_{\left(j\right)})}+{\mathit{b}}_{(k)}$$

(23)

Comparing equation (23) with equations (14) and (15) and pairing the corresponding terms, one can define as

$${\mathit{t}}_{(k)(j)}=\frac{1}{V_{(j)}}\frac{\partial W_{\left(k\right)}}{\partial ({\mathit{u}}_{\left(j\right)}-{\mathit{u}}_{\left(k\right)})}$$

(24a)

and

$${\mathit{t}}_{(j)(k)}=\frac{1}{V_{(k)}}\frac{\partial W_{\left(j\right)}}{\partial ({\mathit{u}}_{\left(k\right)}-{\mathit{u}}_{\left(j\right)})}$$

(24b)

which can also be explicitly written as

$${\mathit{t}}_{(k)(j)}=\left\{\begin{array}{c}{t}_1^{(k)(j)}\hfill \\ \multicolumn{1}{c}{t_2^{(k)(j)}}\\ \multicolumn{1}{c}{t_z^{(k)(j)}}\end{array}\right\}=\frac{1}{V_{(j)}}\left\{\begin{array}{c}\frac{\partial W_{(k)}}{\partial (u_1^{\left(j\right)}-u_1^{\left(k\right)})}\hfill \\ \multicolumn{1}{c}{\frac{\partial W_{(k)}}{\partial (u_2^{\left(j\right)}-u_2^{\left(k\right)})}}\\ \multicolumn{1}{c}{\frac{\partial W_{(k)}}{\partial (w_{\left(j\right)}-w_{\left(k\right)})}}\end{array}\right\}$$

(25a)

and

$${\mathit{t}}_{(j)(k)}=\left\{\begin{array}{c}{t}_1^{(j)(k)}\hfill \\ \multicolumn{1}{c}{t_2^{(j)(k)}}\\ \multicolumn{1}{c}{t_z^{(j)(k)}}\end{array}\right\}=\frac{1}{V_{(k)}}\left\{\begin{array}{c}\frac{\partial W_{(j)}}{\partial (u_1^{\left(k\right)}-u_1^{\left(j\right)})}\hfill \\ \multicolumn{1}{c}{\frac{\partial W_{\left(j\right)}}{\partial (u_2^{\left(k\right)}-u_2^{\left(j\right)})}}\\ \multicolumn{1}{c}{\frac{\partial W_{\left(j\right)}}{\partial (w_{\left(k\right)}-w_{\left(j\right)})}}\end{array}\right\}$$

(25b)

where $u_1$ , $u_2$ and $w$ are the components of the displacement vector.

The strain energy density function given in equation (12) can be written in PD form for the material point $k$ . This can be achieved by transforming all differential terms into the equivalent form of integration. As derived in Appendix A, the partial differential terms in equation (12) for the material point k can be expressed in PD form as

$$u_{i,l}^{\left(k\right)}{u}_{i,l}^{\left(k\right)}+u_{i,l}^{\left(k\right)}{u}_{l,i}^{\left(k\right)}+u_{i,i}^{\left(k\right)}{u}_{l,l}^{\left(k\right)}=\frac{12}{\pi {\delta }^{3}h}{\sum }_{i^k=1}^N\frac{(u_i^{(i^k)}-u_i^{(k)})(u_l^{(i^k)}-u_l^{(k)})}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}$$

(26a)

$$u_{i,i}^{\left(k\right)}{u}_{l,l}^{\left(k\right)}={\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2{\sum }_{i^k=1}^N\frac{u_i^{(i^k)}-u_i^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{j=1}^N\frac{u_l^{(i^k)}-u_l^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}$$

(26b)

$$w_{\mathrm{ii}}^{\left(k\right)}{w}_{\mathrm{ll}}^{\left(k\right)}={\left(\frac{4}{\pi {\delta }^{2}h}\right)}^2{\left({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\right)}^2$$

(26c)

$$\begin{array}{c}{w}_{,\mathrm{il}}^{\left(k\right)}{w}_{,\mathrm{il}}^{\left(k\right)}-w_{\mathrm{ii}}^{\left(k\right)}{w}_{\mathrm{ll}}^{\left(k\right)}=2{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\hfill \\ [8{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{j=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-5{\left({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\right)}^2]\hfill \end{array}$$

(26d)

$$u_{i,i}^{\left(k\right)}{w}_{\mathrm{ll}}^{\left(k\right)}=2{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2{\sum }_{i^k=1}^N\frac{u_i^{(i^k)}-u_i^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}$$

(26e)

$$\begin{array}{c}\frac{1}{2}(u_{i,j}^{(k)}+u_{j,i}^{(k)})w_{,\mathrm{ij}}^{(k)}-u_{i,i}^{(k)}{w}_{,\mathrm{jj}}^{(k)}={\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\hfill \\ (2{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}\hfill \\ +2{\sum }_{i^k=1}^N\frac{u_l^{\left(i^k\right)}-u_l^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}\hfill \\ -3{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)})\hfill \end{array}$$

(26f)

where $n_1=\cos \phi ,n_2=\sin \phi$ with $\phi$ being the bond angle with respect to x-axis. Substituting the relationships given in equations (26a–f) in equation (12), the strain energy density components for the material point $k$ can be written in PD form as

$$\begin{array}{c}{W}_I^{(k)}=\frac{3}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{(1-\nu {\left(z\right)}^2)}\hfill \\ [(1-\nu \left(z\right))\left(\frac{1}{\pi {\delta }^{3}h}{\sum }_{i^k=1}^N\frac{(u_i^{(i^k)}-u_i^{(k)})(u_l^{(i^k)}-u_l^{(k)})}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}\right)\hfill \\ +(\nu \left(z\right)-\frac{1}{3}){\left(\frac{1}{\pi {\delta }^{2}h}\right)}^2{\sum }_{i^k=1}^N\frac{u_i^{(i^k)}-u_i^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{u_l^{(i^k)}-u_l^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}]\mathrm{dz}\hfill \end{array}$$

(27a)

$$\begin{array}{c}{W}_{\mathrm{II}}^{(k)}={\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z^2}{(1-\nu {\left(z\right)}^2)}\hfill \\ [\begin{array}{c}2{\left({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\right)}^2+(1-\nu \left(z\right))\hfill \end{array}\hfill \\ [8{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-5{\left({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\right)}^2]]\mathrm{dz}\hfill \end{array}$$

(27b)

$$\begin{array}{l}{W}_{III}^{(k)}=-{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z}{(1-\nu {\left(z\right)}^2)}\\ [2{\sum }_{i^k=1}^N\frac{u_i^{(i^k)}-u_i^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{(i^k)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\\ +(1-\nu \left(z\right)(2{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{(i^k)}\\ +2{\sum }_{i^k=1}^N\frac{u_l^{\left(i^k\right)}-u_l^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{(i^k)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}\\ -3{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}V(i^k){\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)})]dz\end{array}$$

(27c)

By changing the indices in equations (27a–c), the strain energy density of the material point j, which is a family member of the material point k, can be written as

$$\begin{array}{c}{W}_I^{(j)}=\frac{3}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{(1-\nu {\left(z\right)}^2)}\hfill \\ [(1-\nu \left(z\right))\left(\frac{1}{\pi {\delta }^{3}h}{\sum }_{i^j=1}^N\frac{(u_i^{(i^j)}-u_i^{(j)})(u_l^{(i^j)}-u_l^{(j)})}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{n}_l^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)}\right)\hfill \\ +(\nu \left(z\right)-\frac{1}{3}){\left(\frac{1}{\pi {\delta }^{2}h}\right)}^2{\sum }_{i^j=1}^N\frac{u_i^{(i^j)}-u_i^{(j)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)}{\sum }_{i^j=1}^N\frac{u_l^{(i^j)}-u_l^{(j)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_l^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)}]\mathrm{dz}\hfill \end{array}$$

(28a)

$$\begin{array}{c}{W}_{\mathrm{II}}^{(k)}=\hfill \\ {\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z^2}{(1-\nu {\left(z\right)}^2)}\hfill \\ \{\begin{array}{c}2{\left({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\right)}^2+(1-\nu \left(z\right))\hfill \end{array}\hfill \\ \left[\begin{array}{c}8{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-5{\left({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\right)}^2\hfill \end{array}\right]\}\mathrm{dz}\hfill \end{array}$$

(28b)

$$\begin{array}{l}{W}_{III}^{(k)}=-{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E\left(z\right)·z}{(1-\nu {\left(z\right)}^2)}\\ [2{\sum }_{i^k=1}^N\frac{u_i^{(i^k)}-u_i^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\\ +(1-\nu (z))(2{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}\\ +2{\sum }_{i^k=1}^N\frac{u_l^{\left(i^k\right)}-u_l^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{(i^k)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}\\ -3{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)})]\mathrm{dz}\end{array}$$

(28c)

After substituting equations (27) and (28) into equations (11) and (23) and performing some algebraic manipulations, the PD equation of motion can be written in terms of displacements only as

$$\begin{array}{c}\multicolumn{1}{c}{\rho {\stackrel{··}{u}}_m^{(k)}=\frac{12}{\pi {\delta }^{3}h}\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{1+\nu \left(z\right)}\mathrm{dz}{\sum }_{j=1}^N\frac{u_i^{(j)}-u_i^{(k)}}{{\xi }_{\left(j\right)\left(k\right)}}{n}_i^{\left(j\right)\left(k\right)}{n}_m^{\left(j\right)\left(k\right)}{V}_{(j)}}\\ +{\left(\frac{1}{\pi {\delta }^{2}h}\right)}^2\frac{6}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)(\nu \left(z\right)-\frac{1}{3})}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\sum }_{j=1}^N\frac{n_m^{\left(j\right)\left(k\right)}}{{\xi }_{\left(j\right)\left(k\right)}}\hfill \\ [{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}+{\sum }_{i^j=1}^N\frac{u_i^{\left(i^j\right)}-u_i^{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)}]V_{\left(j\right)}\hfill \\ -2{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E\left(z\right)·z}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\sum }_{j=1}^N\frac{n_m^{\left(j\right)\left(k\right)}}{{\xi }_{\left(j\right)\left(k\right)}}\hfill \\ [{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}+{\sum }_{i^j=1}^N\frac{w_{\left(i^j\right)}-w_{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)}]V_{\left(j\right)}\hfill \\ -{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E\left(z\right)·z}{(1+\nu \left(z\right))}\mathrm{dz}\hfill \\ \{4{\sum }_{j=1}^N\frac{n_l^{\left(j\right)\left(k\right)}}{{\xi }_{\left(j\right)\left(k\right)}}({\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_m^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}+{\sum }_{i^j=1}^N\frac{w_{\left(i^j\right)}-w_{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{n}_m^{\left(i^j\right)\left(j\right)}{n}_l^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)})V_{\left(j\right)}\hfill \\ -3{\sum }_{j=1}^N\frac{n_m^{\left(j\right)\left(k\right)}}{{\xi }_{\left(j\right)\left(k\right)}}({\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}+{\sum }_{i^j=1}^N\frac{w_{\left(i^j\right)}-w_{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)})V_{\left(j\right)}\}+b_m^{(k)}\hfill \\ \mathrm{with}m=1,2\hfill \end{array}$$

(29)

and

$$\begin{array}{l}\rho {\stackrel{··}{w}}_{(k)}\\ =4{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z^2}{(1-\nu {\left(z\right)}^2)}dz{\sum }_j^N=1\frac{1}{{\xi }_{\left(j\right)\left(k\right)}^2}\{{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{w_{(i^j)}-w_{(j)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)}\\ +\frac{1-\nu \left(z\right)}{2}[8n_i^{\left(j\right)(k)}{n}_l^{\left(j\right)\left(k\right)}({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{w_{(i^j)}-w_{(j)}}{{\xi }_{(i^j)\left(j\right)}^2}{n}_i^{\left(i^j\right)\left(j\right)}{n}_l^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)})\\ -5({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{w_{(i^j)}-w_{(j)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)})]\}{V}_{(j)}\\ -{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E\left(z\right)·z}{\left(1-\nu {\left(z\right)}^2\right)}dz{\sum }_j^N=1\frac{1}{{\xi }_{\left(j\right)\left(k\right)}^2}\\ \{2({\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{u_i^{\left(i^j\right)}-u_i^{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)})\\ +\left(1-\nu \left(z\right)\right)[4n_i^{\left(j\right)\left(k\right)}{n}_l^{\left(j\right)(k)}({\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{u_i^{\left(i^j\right)}-u_i^{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_l^{i^j\left(j\right)}{V}_{\left(i^j\right)})\\ -3({\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{u_i^{\left(i^j\right)}-u_i^{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)})]\}{V}_{\left(j\right)}+b_z^{(k)}\end{array}$$

(30)

In the particular case of our study, when the plate is subjected to a constraint condition on all edges, i.e. no free edges, the strain energy density function of a particular material point of the plate can be considerably simplified, as explained in Appendix B, and takes the following form as

$$W=\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{2(1-\nu {\left(z\right)}^2)}\left\{\begin{array}{c}\frac{1-\nu \left(z\right)}{2}(u_{i,j}{u}_{i,j}+u_{i,j}{u}_{j,i}+u_{i,i}{u}_{j,j})+\frac{3\nu \left(z\right)-1}{2}{u}_{i,i}{u}_{j,j}\hfill \\ \multicolumn{1}{c}{+z^2·w_{,\mathrm{ii}}{w}_{,\mathrm{jj}}}\\ \multicolumn{1}{c}{-2z·u_{i,i}{w}_{,\mathrm{jj}}}\end{array}\right\}\mathrm{dz}$$

(31)

with

$$W_I=\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{2(1-\nu {\left(z\right)}^2)}[\frac{1-\nu \left(z\right)}{2}(u_{i,j}{u}_{i,j}+u_{i,j}{u}_{j,i}+u_{i,i}{u}_{j,j})+\frac{3\nu \left(z\right)-1}{2}{u}_{i,i}{u}_{j,j}]\mathrm{dz}$$

(32a)

$$W_{\mathrm{II}}=\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z^2}{2(1-\nu {\left(z\right)}^2)}\mathrm{dz}·w_{,\mathrm{ii}}{w}_{,\mathrm{jj}}$$

(32b)

$$W_{\mathrm{III}}=-\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z}{1-\nu {(z)}^2}\mathrm{dz}·u_{i,i}{w}_{,\mathrm{jj}}$$

(32c)

Corresponding PD form for the material point $k$ can be written as

$$\begin{array}{c}{W}_I^{(k)}=\frac{3}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{(1-\nu {\left(z\right)}^2)}\hfill \\ [(1-\nu \left(z\right))\left(\frac{1}{\pi {\delta }^{3}h}{\sum }_{i^k=1}^N\frac{(u_i^{(i^k)}-u_i^{(k)})(u_l^{(i^k)}-u_l^{(k)})}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}\right)\hfill \\ +(\nu \left(z\right)-\frac{1}{3}){\left(\frac{1}{\pi {\delta }^{2}h}\right)}^2{\sum }_{i^k=1}^N\frac{u_i^{(i^k)}-u_i^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{u_l^{(i^k)}-u_l^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_l^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}]\mathrm{dz}\hfill \end{array}$$

(33a)

$$W_{\mathrm{II}}^{(k)}=2{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z^2}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\left({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}\right)}^2$$

(33b)

$$W_{\mathrm{III}}^{(k)}=-2{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)·z}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\sum }_{i^k=1}^N\frac{u_i^{(i^k)}-u_i^{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}{\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}$$

(33c)

which can also be written concise form as

$$\begin{array}{c}\rho {\stackrel{··}{u}}_m^{(k)}=\hfill \\ \frac{12}{\pi {\delta }^{3}h}\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)}{1+\nu \left(z\right)}\mathrm{dz}{\sum }_{j=1}^N\frac{u_i^{(j)}-u_i^{(k)}}{{\xi }_{\left(j\right)\left(k\right)}}{n}_i^{\left(j\right)\left(k\right)}{n}_m^{\left(j\right)\left(k\right)}{V}_{(j)}\hfill \\ +{\left(\frac{1}{\pi {\delta }^{2}h}\right)}^2\frac{6}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z)(\nu \left(z\right)-\frac{1}{3})}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\sum }_{j=1}^N\frac{n_m^{\left(j\right)\left(k\right)}}{{\xi }_{\left(j\right)\left(k\right)}}\hfill \\ [{\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}+{\sum }_{i^j=1}^N\frac{u_i^{\left(i^j\right)}-u_i^{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)}]V_{\left(j\right)}\hfill \\ -2{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E\left(z\right)·z}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\sum }_{j=1}^N\frac{n_m^{\left(j\right)\left(k\right)}}{{\xi }_{\left(j\right)\left(k\right)}}\hfill \\ [{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}+{\sum }_{i^j=1}^N\frac{w_{\left(i^j\right)}-w_{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)}]V_{\left(j\right)}+b_m^{(k)}\hfill \end{array}$$

(34a)

$$\begin{array}{c}\rho {\stackrel{··}{w}}_{(k)}=\hfill \\ 4{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E(z).z^2}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\sum }_{j=1}^N\frac{1}{{\xi }_{\left(j\right)\left(k\right)}^2}\hfill \\ ({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{w_{(i^j)}-w_{(j)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)})V_{(j)}\hfill \\ -{\left(\frac{2}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\frac{E\left(z\right).z}{(1-\nu {\left(z\right)}^2)}\mathrm{dz}{\sum }_{j=1}^N\frac{2}{{\xi }_{\left(j\right)\left(k\right)}^2}\hfill \\ ({\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{u_i^{\left(i^j\right)}-u_i^{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)})V_{\left(j\right)}+b_z^{(k)}\hfill \end{array}$$

(34b)

Note that if the Poisson’s ratio is considered as a constant value, $\nu (z)=\frac{1}{3}$ , then the PD equations will further simplify as

$$\begin{array}{c}\rho {\stackrel{··}{u}}_m^{(k)}=\frac{9}{\pi {\delta }^{3}h}\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}E(z)\mathrm{dz}{\sum }_{j=1}^N\frac{u_i^{(j)}-u_i^{(k)}}{{\xi }_{\left(j\right)\left(k\right)}}{n}_i^{\left(j\right)\left(k\right)}{n}_m^{\left(j\right)\left(k\right)}{V}_{(j)}-{\left(\frac{3}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}E\left(z\right)·\mathrm{zdz}{\sum }_{j=1}^N\frac{n_m^{\left(j\right)\left(k\right)}}{{\xi }_{\left(j\right)\left(k\right)}}\hfill \\ [{\sum }_{i^k=1}^N\frac{w_{\left(i^k\right)}-w_{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}+{\sum }_{i^j=1}^N\frac{w_{\left(i^j\right)}-w_{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)}]V_{\left(j\right)}+b_m^{(k)}\hfill \end{array}$$

(35a)

$$\begin{array}{c}\rho {\stackrel{··}{w}}_{(k)}=\hfill \\ 2{\left(\frac{3}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}E(z)·z^2\mathrm{dz}{\sum }_{j=1}^N\frac{1}{{\xi }_{\left(j\right)\left(k\right)}^2}\hfill \\ ({\sum }_{i^k=1}^N\frac{w_{(i^k)}-w_{(k)}}{{\xi }_{\left(i^k\right)\left(k\right)}^2}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{w_{(i^j)}-w_{(j)}}{{\xi }_{\left(i^j\right)\left(j\right)}^2}{V}_{\left(i^j\right)})V_{(j)}-{\left(\frac{3}{\pi {\delta }^{2}h}\right)}^2\frac{1}{h}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}E\left(z\right)·\hfill \\ \mathrm{zdz}{\sum }_{j=1}^N\frac{1}{{\xi }_{\left(j\right)\left(k\right)}^2}({\sum }_{i^k=1}^N\frac{u_i^{\left(i^k\right)}-u_i^{\left(k\right)}}{{\xi }_{\left(i^k\right)\left(k\right)}}{n}_i^{\left(i^k\right)\left(k\right)}{V}_{\left(i^k\right)}-{\sum }_{i^j=1}^N\frac{u_i^{\left(i^j\right)}-u_i^{\left(j\right)}}{{\xi }_{\left(i^j\right)\left(j\right)}}{n}_i^{\left(i^j\right)\left(j\right)}{V}_{\left(i^j\right)})V_{\left(j\right)}+b_z^{(k)}\hfill \end{array}$$

(35b)

4. Numerical results

To demonstrate the validity of the presented PD formulation for functionally graded Kirchhoff plates, the PD solutions are compared with the corresponding finite element (FE) analysis results. Here, the material properties are chosen such that Young’s modulus, $E(z)$ , is assumed to be varied linearly through the thickness direction and Poisson’s ratio, $\nu (z)$ , remains a constant as

$$E\left(z\right)=(E_t-E_b)\frac{z}{h}+\frac{1}{2}(E_t+E_b)\left(\mathrm{GPa}\right)$$

(36a)

and

$$\nu (z)=\nu$$

(36b)

where $E_t$ and $E_b$ denote the Young’s modulus of the top and bottom surfaces of the plate, and $h$ represents the total thickness of the plate.

4.1. Clamped plate subjected to transverse loading

A clamped functionally graded plate with a length and width of $L=W=1\mathrm{m}$ and a thickness of $h=0.02\mathrm{m}$ is considered, as shown in Figure 2. The Poisson’s ratio of the plate is $\nu =1/3$ and the Young’s modulus of the top and bottom surfaces are $E_t=200\mathrm{GPa}$ and $E_b=100\mathrm{GPa}$ , respectively. The model is discretised into one single row of material points along the thickness direction and the distance between material points is $\mathrm{\Delta }x=1/101\mathrm{m}$ . A fictitious region is introduced outside the edges as the external boundaries with a width of $2\delta$ . The plate is subjected to a distributed transverse load of $p_z=101\mathrm{N}/\mathrm{m}$ through the y-centre line, respectively. The line load is converted to a body load of $b=\frac{p_{z}W}{(\frac{W}{\mathrm{\Delta }x})\mathrm{\Delta }V}=510050\mathrm{N}/{\mathrm{m}}^3$ and it is distributed to one column of material volumes through the central line.

Figure 2.

Clamped plate subjected to transverse loading.

The FE model of the plate is created by using the SHELL181 element in ANSYS with dimensions of $1\mathrm{m}\times 1\mathrm{m}\times 0.02\mathrm{m}$ . To model the functionally graded plate, the model is divided to 50 layers with varying homogeneous materials properties through the thickness. The Young’s modulus varies linearly over the thickness from the first layer $E_1=101\mathrm{GPa}$ to the last layer $E_{50}=199\mathrm{GPa}$ , as shown in Figure 3. The model is meshed with 0.01 m element size.

Figure 3.

Variation of the Young’s modulus in thickness direction for the FE model.

The PD and FE transverse displacement contours are compared in Figure 4. They yield similar displacement variations. The maximum difference between the PD and FE results is less than 0.5%. Moreover, the transverse and in-plane displacement components along the central y-axis are compared in Figure 5 and very good agreement is obtained between the two approaches. These results verify the accuracy of the current PD formulation for a functionally graded Kirchhoff plate theory under clamped boundary conditions.

Figure 4.

Variation of transverse displacements (a) PD, (b) FEM.

Figure 5.

Variation of (a) transverse, w and (b) in-plane, u₁, displacements along the central y-axis.

4.2 Simply supported plate subjected to transverse loading

A simply supported plate (see Figure 6) has the same geometrical and material properties as in the clamped plate case. Again, it is discretised with a single row of material points along the thickness direction and the discretisation size is $\mathrm{\Delta }x=\frac{1}{101}m$ . A fictitious region is created outside the region of boundaries and its width is equal to two times the size of the horizon, $\delta$ . The plate is subjected to a distributed transverse line load of $p_z=101N/m$ through the y-central line. It is imposed to central column of material points with a body load of $b=\frac{p_{z}W}{(\frac{W}{\mathrm{\Delta }x})\mathrm{\Delta }V}=510050N/m^3$ .

Figure 6.

Simply supported plate subjected to transverse loading.

The transverse displacement components of FE and PD theory show very close variations, as shown in Figure 7. The maximum difference between the PD and FE results is less than 0.5%. Furthermore, the transverse and in-plane displacement variations along the central y-axis are on top of each other for the FE and PD results, as shown in Figure 8. This confirms the current PD formulation of functionally graded plate theory under simply supported boundary conditions.

Figure 7.

Variation of transverse displacements (a) PD, (b) FEM.

Figure 8.

Variation of (a) transverse, w and (b) in-plane, u₁, displacements along the central y-axis.

4.3 Simply supported plate subjected to inclined loading

This problem case is similar to the previous case except Poisson’s ratio of 0.2 and an inclined load of $p_x=p_z=101N/m$ through the y-central line are considered, as shown in Figure 9.

Figure 9.

Simply supported plate subjected to inclined loading.

As depicted in Figure 10, the transverse displacement components of FE and PD theory agree very well with each other. Moreover, the transverse and in-plane displacements variations along the central y-axis are also in very good agreement, as depicted in Figure 11.

Figure 10.

Variation of transverse displacements (a) PD, (b) FEM.

Figure 11.

Variation of (a) transverse, w and (b) in-plane, u₁, displacements along the central y-axis.

5. Conclusions

In this study, a state-based peridynamic model was presented for functionally graded Kirchhoff plates. Equations of motion of the new formulation were obtained by using Euler–Lagrange equations and Taylor’s expansion. The formulation was verified by considering several benchmark problems including a clamped plate subjected to transverse loading and simply supported plate subjected to transverse loading and inclined loading. The material properties were chosen such that Young’s modulus is assumed to be varied linearly through the thickness direction and Poisson’s ratio is constant. Peridynamic results were compared against FE analysis results and a very good agreement is obtained between the two approaches.