Nonlocal dynamic Kirchhoff plate formulation based on nonlocal operator method

In this study, we propose a nonlocal operator method (NOM) for the dynamic analysis of (thin) Kirchhoff plates. The nonlocal Hessian operator is derived based on a second-order Taylor series expansion. The NOM does not require any shape functions and associated derivatives as ’classical’ approaches such as FEM, drastically facilitating the implementation. Furthermore, NOM is higher order continuous, which is exploited for thin plate analysis that requires C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document} continuity. The nonlocal dynamic governing formulation and operator energy functional for Kirchhoff plates are derived from a variational principle. The Verlet-velocity algorithm is used for the time discretization. After confirming the accuracy of the nonlocal Hessian operator, several numerical examples are simulated by the nonlocal dynamic Kirchhoff plate formulation.


Introduction
As a common engineering structure, plate/shell is widely used in civil engineering, aerospace and other fields. The mechanical analysis of rectangular thin plates has always been one of the research focuses of scholars and engineers. The governing equation of the Kirchhoff plate bending problem is a fourth-order partial differential equation whose deflection is an independent variable. The numerical method for developing this problem has a wide scientific significance for solving plate/shell problems.
The analysis of Kirchhoff plate bending problems poses challenges to 'classical' finite element formulations [1][2][3][4][5], which are only C 0 continuous. However, the Kirchhoff plate problem is a fourth-order partial differential equation which requires a C 1 formulation if weak form based methods such as FEM are employed. An efficient alternative are so-called meshless methods as many of them are higher order continuous. The approximation function of the meshless method [6][7][8][9][10][11] represented by the Element-free Galerkin method EFG [12][13][14][15] is highly smooth, that is, and the high-order continuity is satisfied, meantime the meshless method is easy to form high-order. The approximate function has significant advantages to the numerical solutions for higherorder partial-differential equations.
Rabczuk [16] devised a meshfree method for thin shell analysis for finite strains and arbitrary evolving cracks exploiting the higher order continuity of the EFG shape functions and avoiding any rotational degrees of freedom. Mohammed et al. [17,18] presented a meshless method to analyze the mechanical response of elastic thin plates. However, since meshless shape functions are commonly rational functions, more integration points are needed to evaluate the weak form. For example, Brebbia [6] employed 6 × 6 quadrature points, which significantly reduces the computational efficiency. An interesting alternative to meshless methods is isogeometric analysis (IGA) [19,20], which also fulfills the higher order continuity requirement needed for thin plate analysis. This method takes advantage of NURBS/B-Spline basis functions, which are commonly used in computeraided-design (CAD) to describe the geometries. IGA has been successfully applied to the analysis of plates and shells for instance in [21][22][23][24][25]. As CAD geometries are surface representations, they are particularly suited for plates and shells.

Outline of nonlocal operator method
We consider a Kirchhoff plate occupying a domain as illustrated in Fig. 1. Let us denote the spatial coordinates with i , ij ∶= j − i is relative position vector from i to j ; w i ∶= w( i , t) and w j ∶= w( j , t) are the displacement value for i and j , respectively; the relative displacement area for the spatial vector ij is w ij ∶= w j − w i . The support and the dual-support are two fundamental principles in NOM. The domain, in which every spatial point j forms the spatial vector ij from i to j is called the support S i of point i , see Fig. 1b. We can write S = { 2 , 3 , 5 , 6 } . The dual-support of i is defined as a union of points which supports include , i.e.
The dual-support of point j forms the dual-vector i and is denoted as S � = { 3 , 6 } ; ′ ij is the S j relative position space vector. NOM requires fundamental nonlocal operators that replace the local operators in calculus. Thus, the functional designed to construct a residual and tangent stiffness matrix is formulated in terms of the nonlocal differential operator. The higher order nonlocal operator ̃w i for the scalar field w in support S i can be expressed as [34] (1) where ( ij ) represents the weight function, ij represents relative position vector. ij represent the list of polynomials. For example, the polynomials and the higher order nonlocal operator in 2D with maximal second-order derivatives are Using nodal integration, the higher order nonlocal operator and it's variation for a vector field w in discrete forms are We use a penalty energy functional to obtain the linear field of the scalar field to eliminate zero-energy modes. The higher order operator energy functional for a scalar field w at a point i is defined as is the normalization coefficient and w is the penalty coefficient. The nonuniform aspect of the deformation is defined as T ij̃w i − w ij , and the stability of the NOM is greatly enhanced while T ij̃w i − w ij is enforced explicitly.

Derivation of nonlocal Hessian operator for Kirchhoff plate
In this part, we derive the second-order nonlocal Hessian operator and its variation. The second derivative and its variation in 1D is given as To devise the nonlocal Hessian operator in 2D, let us define first the vector ij = (x ij , y ij ) T . It can be shown that the second-order shape tensor for point i is computed by [30] (2) x 2 ij x ij y ij x ij y ij y 2 ij dV j and the third-order shape tensor for point i by which finally leads to The fourth-order shape tensor for point i is computed by The second-order Taylor series extension for a scalar field w is given as so that we obtain which finally results in the following Hessian nonlocal operator and (9) The weighted tensor 1 2 ∇ T ∇w ∶ 4i can be simplified to Note that the rank of 4i is 3 and the 2D nonlocal Hessian operator has only three independent variables where 4i as the pseudo-inverse of 4i in this study.
For the scalar field w, the 2D nonlocal Hessian operator for point i can be written in matrix form as 2i and the pseudo-inverse of 4i be explicitly written as To facilitate the calculation −1 4i , we convert 4i to a 3 × 3 matrix yielding 2i ij is a matrix 2×2 with the terms Q 12 = Q 21 . Remove the repeated terms in matrix , and the remain terms Q 11 , Q 12 , Q 22 in matrix can be reconstituted a vector 3×1 = (Q 11 , Q 12 , Q 22 ) T . Then Eq.16 can be rewritten as can be obtained by reconstituting terms in Eq.21 and can be expressed as where The explicit form of the nonlocal Hessian operator in 2D can finally be expressed by Note that e i11 , e i12 , e i22 are calculated for each neighbor in the support domain. As a result, the variation of nonlocal Hessian operator can be obtained in explicit form as 3 Derivation of nonlocal dynamic Kirchhoff plate formulation

Classical elastic plate theory
Kirchhoff plate theory assumes that the normal stress in the thickness direction can be ignored and the normal of the midplane of the plate remains normal after deformation. Hence, all stresses and strains can be expressed by the deflection w of the midplane of the plate. Considering the plate element shown in Fig. 2, the in-plane displacements u and v can, therefore, be expressed in terms of the first derivatives of w, i.e.
The strain resultant can be obtained by where x and y indicates the curvature of the midplane of the plate in the x-and y-direction while xy refers to the torsion, respectively. The stress resultants of the Kirchhoff plate are given by where M x and M y are the bending moment per unit length around the y-and negative x-axes, respectively, while M xy (= M yx ) is the torque per unit length. With a linear stress distribution in the z-direction and assuming a thickness of t, the stresses can be computed by The Cauchy stress tensor can also be expressed in terms of the linear strain tensor assuming Hooke's law: y).
with Finally, the constitutive model can be formulated in terms of the stress and strain resultants by where plate is the constitutive matrix defined as is the Kirchhoff plate's bending stiffness, such that Eq.(30) can be rewritten as

Nonlocal dynamic Kirchhoff plate formulation
The total Lagrange energy functional for the Kirchhoff plate can be expressed as with ẇ = w t ; is the density of the plate and q z a distributed load in the z-direction. Replacing the local Hessian ∇ T ∇w with the nonlocal Hessian ∇ T∇ w in Eq.34, we obtain 12 21 22 . where the boundary condition w(t 1 ) = 0, w(t 2 ) = 0 is considered in the above derivation. According to Hamilton's principle, for any w i , the first variation of the functional ϝ should be zero, which leads to The nonlocal form is correlated to the local form by According to Eq.23, we devise the explicit form of ∇ T∇ ∶ i As Eq. 37 suffers from zero-energy modes, we introduce the so-called nonlocal operator energy functional, which is described in the next section.

Operator energy functional
For the Kirchhoff plate, the maximal order of partial derivatives in Eq.37 is two, hence we select the second order of nonlocal operators in Eq.5. The operator energy functional for second order nonlocal operators of a scalar field w for point i can be expressed as e j21 e j22 dV j . For the scalar field w, the internal force due to the operator energy functional is given by indicates the zero-energy internal force. Finally, the correspondence between local and nonlocal formulation and the operator functional enhanced governing equation of Kirchhoff plate can be expressed as Similarly, parallel to the x− axis (as y = b simply supported boundary), the boundary conditions is

Numerical implementation
The domain Ω is decomposed into N points occupying a volume ΔV i : For each point, the support is denoted by where j 1 , ..., j k , ..., j n i are the global indices of the neighbors of point i and n i represents the number of neighbors in support domain S i . The out-of-plane transversal force between points are computed by

Kirchhoff plate boundary conditions
Let us consider the boundary conditions shown in Fig. 3, which can be classified into: 1. Clamped boundary conditions, where the deflection and the slope of the mid-plane is zero. The positioning shift w and the section rotation are both zero. The boundary conditions are in the direction parallel to the y-axis (as x = 0 clamped boundary): Parallel to the x−axis (as y = 0 clamped boundary), the Kirchhoff plate boundary conditions is 2. Simply supported Kirchhoff plate boundary conditions where the plate is free to rotate about a line but prevented from deflecting. The positioning shift w and moment M n value is zero: parallel to the y−axis (as x = a simply supported boundary), the boundary conditions is and Eq.46 can be written as (45) Finally, we get the internal force P ij between points as So that in discrete form, Newton's equation of motion is expressed as where t is the time, w i (t) = w 1 (t), … , w N (t) is the ensemble of the position vector of N points, F i is the external force vector and P i denotes the internal force vector; i refers to the mass of the point. In this paper, the velocity and displacement is updated via the Verlet-Velocity scheme [35]: For reasons of stability, we applied a damping term to each point F s i representing the damping force for each point, ẇ i represents the velocity (vector) of the point, and c is a damping coefficient.
The main implementation process of higher order explicit NOM for the dynamic analysis of Kirchhoff plate can be summarized as follows:

Numerical examples
The proposed nonlocal dynamic Kirchhoff plate formulation is implemented in Wolfram Mathematica. After verifying the accuracy of the nonlocal Hessian operator, several benchmark problems are studied and compared to results obtained by ABAQUS using the S4R plate/shell element [36].

Verification of nonlocal Hessian operator
Let us consider a simply supported square Kirchhoff plate with a width of a 0 = 10 m and thickness t=0.01m. The plate is subjected to a uniform pressure of q z = 100N∕m 2 . Young's modulus and Poisson's ratio are E=30GPa and =0.3, respectively. The number of neighbors of each point is set to n = 24 . To test the accuracy of the nonlocal Hessian operator, we assume Υ = 1, 3 (see Eq.58). The analytical solution of this problem is given in [37], i.e.
with Υ = Υ 2 . We employ the quintic spline function as weight function To accurately represent the relative error of the operator, we consider points at y = 0 . Figs. 5, 6 show the deflection curve of the numerical simulation (

Nonlocal dynamic Kirchhoff plate formulation with simply supported boundary condition
We now focus on a simply supported square Kirchhoff plate with a width of a 0 = 1 m and thickness t=0.01m. Young's modulus and Poisson's ratio are E=200GPa and =0.3, respectively. A uniform pressure of q z = 100N∕m 2 is applied to the plate. The number of neighbors in the domain of influence of each point is set to n = 24 . The distance between points is selected as Δx=0.01 m leading to 10201 points. For the ABAQUS model-as a comparison-we discretized the plate into 100 × 100 elements using the same input parameters.
Simply supported boundary conditions are assumed: Contour plots of the deflection can be found in Figs. 7 and 8, respectively.

Nonlocal dynamic Kirchhoff plate formulation with clamped boundary condition
Next, a plate using clamped boundary conditions is studied. The geometry and the parameters of the previous section are adopted. The clamped boundary conditions are given by Contour plots of the deflection are illustrated in Figs. 9 and 10, respectively.

Transient test of nonlocal dynamic Kirchhoff plate formulation with simply supported boundary condition
The last example is a simply supported plate where the damping is omitted. A uniform pressure of q z = 100N∕m 2 is applied to the plate and the neighbors assigned to each point is set to n = 30 . All other parameters are adopted from the previous example. Figure 11 shows the deflection field at different times.

Conclusions
In this paper, a nonlocal dynamic Kirchhoff plate formulation based on a nonlocal operator method is proposed. Therefore, we derived the explicit form of the nonlocal Hessian operator taking advantage of a second order Taylor series expansion. The nonlocal operator energy functional is also derived and the relationship between local and nonlocal formulations is interpreted. In the numerical simulation section, we first verify the accuracy of the nonlocal Hessian operator for the Kirchhoff plate and compare it to analytical solutions. Subsequently, the nonlocal dynamic Kirchhoff plate formulation with different type boundary conditions (clamped and simply supported) is studied and compared to simulations obtained by ABAQUS. In the future, we intend to extend the formulation for nonlinear dynamic fracture exploiting the advantages and flexibility of NOM in this direction.
Acknowledgements The research was supported by the China Scholarship Council. The author would like to express their sincere gratitude to the editor and two anonymous reviewers for their valuable comments, which further enhanced the quality of this manuscript (62) w(x, 0) = w(x, 1) = w(0, y) = w(1, y) = 0 w(0, y) Funding Open Access funding enabled and organized by Projekt DEAL.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.