A Symplectic Method for Analyzing the Nonlocal Modal Behavior of Kirchhoff Plates and Numerical Validation

1. Introduction

In recent years, there has been a rapid expansion in the technical application of MEMS/NEMS nanoplates in engineering fields. From ultra-sensitive mass sensors and high-frequency radio frequency (RF) filters in the industrial sector to single-molecule detection chips and DNA sequencing devices in the biomedical field, there is an increasing shift from traditional macroscopic structures to nanostructures. This shift is driven in part by the high surface-area-to-volume ratio and favorable mechanical properties of nanostructures, which motivates ongoing technological development in related fields. Common nanostructures include nanorods, nanobeams, and nanoplates. It is noteworthy that at the nanoscale, their thermal, electrical, mechanical, physical, and chemical properties often change [1,2,3,4]. Consequently, there is a need to establish specialized research methodologies to address the issues arising from these changes. The following methods may be employed to study the behavior of common nanostructures: experimental analysis [5], molecular dynamics simulations [6], and non-classical continuum mechanics methods [7]. It is evident that, owing to the intricacy of experimental analysis and the time-consuming and costly nature of molecular dynamics simulations, the utilization of non-classical continuum mechanics methods to simulate the mechanical behavior of structures at the nanoscale has become the preferred choice for many researchers.

The more general non-classical continuum mechanics theory was proposed by Eringen and his colleagues [8,9]. In this theory, the stress distribution at a reference point depends on the strain components at that point and the strains at other points in the medium. Furthermore, a kernel function is introduced to simulate the attenuation of long-range interactions between two particles as a function of their distance. The nonlocal elasticity theory under consideration comprises both integral and differential constitutive relations. In comparison with integral constitutive equations, which necessitate consideration of global effects and consequently yield integro-differential equations that demand complex calculations, researchers have historically favored differential constitutive equations for nonlocal studies [10,11,12,13]. However, given that the differential-type constitutive relation is derived under specific conditions, its applicability is limited [14]. Moreover, Eringen's nonlocal theory was originally founded on an integral-type constitutive relation. Consequently, the differential-type is lacking in terms of broader significance and applicability when compared to the integral-type constitutive relation.

Integral constitutive relations have been successfully applied to problems of nonlocal nanobeams and nonlocal Euler beams, resolving issues arising from the use of differential constitutive relations [15,16,17,18]. Compared with one-dimensional structures, nanoplate structures are more complex, and their governing equations are more challenging to handle. Consequently, irrespective of the constitutive relation employed, research on the free vibration of nanoplates is relatively scarce. The following section mentions some representative studies.

Lu et al. proposed Kirchhoff and Mindlin plate element models based on the nonlocal theory [19]. In the subsequent analysis of nonlocal plate vibration, researchers have adopted a wide range of numerical and analytical solution methods. These include the Rayleigh-Ritz method [20,21], the Finite Element Method (FEM) [22,23,24], the Finite Strip Method [25,26], the Galerkin method [27,28], the Differential Quadrature (DQ) method [29,30,31,32], the meshless kp-Ritz method [33], the Chebyshev Collocation Method [34], the Galerkin Strip Distributed Transfer Function Method [35], and the Discrete Singular Convolution Method [36,37].

Numerical methods are of paramount importance in the analysis of the mechanical behavior of nonlocal plates. Inter-belt analysis, a numerical method, was first proposed by Zhong and co-workers [38,39]. Conventional local finite element methods generally consider the connection between elements, or between an element and the external environment, as the displacement at their interface, where the interface is assumed to have no thickness. However, in the context of nonlocal theory, the influence of long-range interactions between elements must be taken into account. In this specific context, the connection zone between two non-adjacent yet mutually influencing elements is no longer a surface but a strip of elements. Since its boundary is no longer a plane but has a finite width, it is termed an inter-belt. Since its proposal, this analysis method has been applied to the nonlocal analysis of carbon nanotubes [40] and subsequently verified for nonlocal rods, Euler beams, Timoshenko beams, and plane-strain problems [41,42,43].

This study applies the inter-belt method to vibration analysis of nonlocal Kirchhoff plates to examine the effects of different nonlocal parameters, two-phase mixing parameters, mode numbers, kernel functions, and geometric parameters on vibration frequencies under nonlocal conditions. In addition, the distinctive advantages of this algorithm in the resolution of integro-differential equations are substantiated.

The subsequent section of this paper will provide a review of the fundamental form of Eringen's nonlocal elasticity theory and discuss the methodology for handling integral constitutive relations. In the subsequent Section 3, the nonlocal finite element solution for Kirchhoff plates is implemented. The fourth section of this text presents numerical results and conducts a comprehensive parametric study. The findings of the present study will be discussed in Section 5, while the conclusions will be presented in Section 6.

2. Nonlocal Elasticity Theory

2.1. Eringen's Nonlocal Elasticity Theory

In accordance with the nonlocal elasticity theory proposed by Eringen [8,9], the stress at a point $x$ in a continuum is contingent not solely on the strain at that particular point, but also on the strains of all other points $x\mathrm{’}$ within the entire body. A constitutive relation in integral form describes this nonlocal effect. In the case of a homogeneous, isotropic linear elastic solid, the constitutive equation can be expressed as

$${\mathit{t}}_{kl}\left(x\right)={\int }_V\alpha \left(\right|x^{\prime }-x|,\kappa ){\mathit{\sigma }}_{kl}\left(x^{\prime }\right)\mathrm{d}v\left(x^{\prime }\right)，$$

(1)

where ${\mathit{t}}_{kl}$ is the nonlocal stress tensor; ${\mathit{\sigma }}_{kl}$ is the classical local macroscopic stress tensor, given by the generalized Hooke's Law.

$${\mathit{\sigma }}_{kl}\left(x^{\prime }\right)=\lambda {\mathit{e}}_{rr}\left(x^{\prime }\right){\delta }_{kl}+2\mu {\mathit{e}}_{kl}\left(x^{\prime }\right)，$$

(2)

where $\lambda$ and $\mu$ are Lame constants, ${\mathit{e}}_{kl}$ is the linear strain tensor, and ${\delta }_{kl}$ is the Kronecker delta.

The equation of motion reads

$${\mathit{t}}_{kl,k}+\rho (f_l-{\ddot{\mathit{u}}}_l)=0.$$

(3)

The kernel function $\alpha \left(\right|x^{\prime }-x|,\kappa )$ describes the attenuation characteristics of nonlocal interactions, where $\kappa =e_0{a}_0/l_0$ is the dimensionless nonlocal parameter, $a_0$ is the internal characteristic length (e.g., lattice constant), $l_0$ is the external characteristic length (e.g., wavelength or crack length), and $e_0$ is a constant determined by material properties.

Although the integral constitutive relation can be transformed into a differential form for specific kernel functions, the integral form is more broadly applicable. It can handle more complex boundary conditions and physical phenomena.

2.2. Inter-Belt Analysis Theory

The inter-belt analysis method is a numerical analysis method developed based on symplectic mathematical theory and computational structural mechanics [38,39]. In traditional finite element theory, the connection between elements is generally regarded as an interface of zero thickness, for which only $C^0$ continuity of displacement is required. However, in the context of nonlocal elasticity theory, the presence of long-range forces extends the interaction between material points beyond the confines of adjacent differential elements, encompassing a finite region termed the influence domain. In this context, the boundary between elements is no longer an idealized geometric surface, but rather a belt-like region with a certain width, i.e., an inter-belt.

In the context of integral nonlocal constitutive relations [41], the inter-belt analysis method has been demonstrated to be an effective solution to the numerical solution problem of integro-differential equations by discretizing the continuous system. In the context of deformation energy involving nonlocal interactions, inter-belt analysis discretizes it into the sum of energies of a series of substructures.

In the specific case of discretizing interactions within the influence domain, the nonlocal global stiffness matrix will no longer be a banded, sparse matrix as in traditional finite elements. Rather, it will be a generalized stiffness matrix that includes coupling terms between non-adjacent nodes. The construction of an element stiffness matrix that incorporates inter-belt effects facilitates the transformation of complex integro-differential dynamic equations into standard linear algebraic eigenvalue problems. This approach enables simulation of the dynamic response of nonlocal structures.

3. Kirchhoff Plate Element and Nonlocal Finite Element Formulation

In the preceding two sections, the fundamental form of Eringen's nonlocal elasticity theory was reviewed, and the rudiments of inter-belt analysis in the discretization of integral nonlocal constitutive relations were introduced. It is on this basis that the present section commences with a concise review of the classical Kirchhoff thin plate element and its local finite element formulation. Subsequently, Eringen's two-phase nonlocal theory is applied to the Kirchhoff plate bending problem. The nonlocal strain energy expression based on the center-point approximation is presented, and the corresponding element stiffness matrix form is derived, providing a theoretical basis for subsequent nonlocal inter-belt analysis, discretization, and numerical examples.

3.1. Classical Kirchhoff Plate Bending Theory and Local Finite Element

The kinetic energy formulation of a classical local Kirchhoff plate is

$$T_e={\displaystyle \frac{1}{2}}{\dot{\mathit{q}}}_e{\left(t\right)}^{\mathrm{T}}\mathit{m}{\dot{\mathit{q}}}_e\left(t\right).$$

(4)

where $\mathit{m}$ is the mass matrix of the plate element, and ${\mathit{q}}_e$ is the element degree-of-freedom vector.

$$\mathit{m}=\rho h{\int }_{-b}^b{\int }_{-a}^a{\mathit{N}}_e{\left(z\right)}^{\mathrm{T}}{\mathit{N}}_e\left(z\right)d\zeta d\eta ,$$

(5)

Assembling the local element stiffness matrices ${\mathit{k}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}$ yields the local global stiffness matrix ${\mathit{K}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}$. Similarly, assembling the local element mass matrices $\mathit{m}$ yields the global mass matrix, $\mathit{M}$. The assembled free-vibration equation reads

$$\mathit{M}\ddot{\mathit{q}}+{\mathit{K}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}\mathit{q}=0$$

(6)

The characteristic equation of the above formula is

$$det\left({\widehat{\mathit{K}}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}-{\omega }^2\widehat{\mathit{M}}\right)=0$$

(7)

where $\widehat{\mathit{K}}$ and $\widehat{\mathit{M}}$ denote the reduced matrices after imposing boundary conditions. The natural frequencies of the local Kirchhoff plates can be obtained by solving the above equation [44]. This study imposes boundary conditions by eliminating rows and columns from the stiffness matrix, and solves the eigenvalue problem using MATLAB.

Under the classical local Kirchhoff plate theory, generalized stress is introduced [45]

$$\mathit{F}={\mathit{D}}_0\mathit{\zeta }={\displaystyle \frac{E{t}^3}{12(1-{\nu }^2)}}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& {\displaystyle \frac{1-\nu }{2}}\end{array}\right]\left[\begin{array}{c}{\displaystyle \frac{\partial {\theta }_y}{\partial x}}\\ -{\displaystyle \frac{\partial {\theta }_x}{\partial y}}\\ -({\displaystyle \frac{\partial {\theta }_y}{\partial y}}+{\displaystyle \frac{\partial {\theta }_x}{\partial x}})\end{array}\right].$$

(8)

The local Kirchhoff plate element stiffness matrix can be obtained through the principle of minimum potential energy as

$${\mathit{k}}_{local}={\int }_{A_E}{\mathit{B}}^{\mathrm{T}}(x,y) {\mathit{D}}_0 \mathit{B}(x,y) \mathrm{d}A$$

(9)

where $\mathit{B}(x,y)\in {\mathbb{R}}^{3\times 12}$ is the geometric matrix of the Kirchhoff plate element. The above theory provides a reference local model for the subsequent introduction of nonlocal effects.

3.2. Derivation of Nonlocal Stiffness for Kirchhoff Plates Based on Eringen's Two-Phase Nonlocal Theory

To characterize the long-range interactions at the nanoscale, the Eringen nonlocal model introduced in Section 2.1 is applied to the Kirchhoff plate bending problem. For a rectangular plate divided into several standard Kirchhoff plate elements, let one non-boundary element be $E0$, with dimensions $2b$ and $2c$ in the $x$ and $y$ directions, respectively, and thickness $t$. Within the framework of the generalized stress-strain relationship, the nonlocal strain energy of element $E0$ is written as

$$U_{nonlocal}={\displaystyle \frac{1}{2}}{\int }_{A_{E0}}{\mathit{\zeta }}^{\mathrm{T}}\left(x,y\right)\left[{\mathit{D}}_0{\int }_{A_{plate}}\alpha \left(s,\kappa \right) \mathit{\zeta }\left(x^{\prime },y^{\prime }\right) \mathrm{d}{A}_{plate}\right]\mathrm{d}{A}_{E0}，$$

(10)

where $A_{plate}$ is the area of the middle surface of the entire plate, $s=\sqrt{(x^{\prime }-x{)}^2+(y^{\prime }-y{)}^2}$ is the two-dimensional Euclidean distance between two integration points (in this paper, it is the geometric center distance of two influencing elements). $\alpha (s,\kappa )$is a nonlocal kernel function that satisfies normalization and monotonic decay conditions, and $D_0$ is the bending stiffness matrix defined in the previous subsection. Discretizing the global integral by elements yields as

$$D_0{\int }_{A_{plate}}\alpha (s,\kappa ) \mathit{\zeta }(x^{\prime },y^{\prime }) \mathrm{d}{A}_{plate}=D_0\sum _{k=0}^n{\int }_{A_{Ek}}{\alpha }_k \mathit{\zeta }(x^{\prime },y^{\prime }) \mathrm{d}{A}_{Ek},$$

(11)

where $A_{Ek}$ is the area of the $k$ -th element. ${\alpha }_k$ is the value of the specific kernel function $\alpha (s,\kappa )$ after substituting $s$ under a fixed $\kappa$, reflecting the strength of long-range interaction between elements. To further simplify, the generalized strain integral over each element is approximated as the volume of a quadrangular prism, with the strain value at the geometric center as the height and the element area as the base.

$${\int }_{A_{Ek}}\mathit{\zeta }(x^{\prime },y^{\prime }) \mathrm{d}{A}_{Ek}\approx {\mathit{\zeta }}_{Ek} A_{Ek},$$

(12)

where ${\mathit{\zeta }}_{Ek}$ is the generalized strain at the geometric center point of element $Ek$. For a standard Kirchhoff rectangular element, it can be expressed through the center point geometric operator $\mathit{T}$,

$${\mathit{\zeta }}_{Ek}=B\left(\mathrm{0,0}\right) {\mathit{q}}_{Ek}=\mathit{T} {\mathit{q}}_{Ek},$$

(13)

where $B\left(\mathrm{0,0}\right)$ is the geometric matrix at the origin of the element's natural coordinates. For a standard element of size $2b\times 2c$, $\mathit{T}$ is explicitly written as

$$\mathit{T}=\mathit{B}\left(\mathrm{0,0}\right)=\left[\begin{array}{cccccccccccc}0& 0& {\displaystyle \frac{1}{4b}}& 0& 0& -{\displaystyle \frac{1}{4b}}& 0& 0& -{\displaystyle \frac{1}{4b}}& 0& 0& {\displaystyle \frac{1}{4b}}\\ 0& -{\displaystyle \frac{1}{4c}}& 0& 0& -{\displaystyle \frac{1}{4c}}& 0& 0& {\displaystyle \frac{1}{4c}}& 0& 0& {\displaystyle \frac{1}{4c}}& 0\\ -{\displaystyle \frac{1}{bc}}& {\displaystyle \frac{1}{4b}}& -{\displaystyle \frac{1}{4c}}& {\displaystyle \frac{1}{bc}}& -{\displaystyle \frac{1}{4b}}& -{\displaystyle \frac{1}{4c}}& -{\displaystyle \frac{1}{bc}}& -{\displaystyle \frac{1}{4b}}& {\displaystyle \frac{1}{4c}}& {\displaystyle \frac{1}{bc}}& {\displaystyle \frac{1}{4b}}& {\displaystyle \frac{1}{4c}}\end{array}\right].$$

(14)

Based on this, it can be concluded as

$${\mathit{D}}_0{\int }_{A_{plate}}\alpha (s,\kappa ) \mathit{\zeta }(x^{\prime },y^{\prime }) \mathrm{d}{A}_{plate}\approx {\mathit{D}}_0\sum _{k=0}^n{\alpha }_k {\mathit{\zeta }}_{Ek} A_{Ek}={\mathit{D}}_0\sum _{k=0}^n{\alpha }_k \mathit{T} {\mathit{q}}_{Ek} A_{Ek}.$$

(15)

In the event of the plate being divided equally into a number of standard rectangular elements of uniform size, the following equality exists for their respective areas,

$$A_{Ek}=A_{E1}=\dots =A_{En}=A_E=4bc.$$

(16)

Substituting the aforementioned results back into the nonlocal strain energy expression for element $E0$.

$$U_{\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}={\displaystyle \frac{1}{2}}{\int }_{A_{E0}}{\mathit{\zeta }}^{\mathrm{T}}\left(x,y\right)\left[{\mathit{D}}_0\sum _{k=0}^n{\alpha }_k \mathit{T} {\mathit{q}}_{Ek}\right]\mathrm{d}{A}_{E0},$$

(17)

$$\mathit{\zeta }\left(x,y\right)=\mathit{B}\left(x,y\right) {\mathit{q}}_{E0}.$$

(18)

Substituting Eqs. (18) and (15) into Eq. (10) yields,

$$U_{\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}=\sum _{k=0}^n{\displaystyle \frac{A_E {\alpha }_k}{2}}{\int }_{A_{E0}}{\mathit{q}}_{E0}^{\mathrm{T}}{B}^{\mathrm{T}}\left(x,y\right) D_0 \mathit{T} {\mathbf{q}}_{Ek} \mathrm{d}{A}_{E0}.$$

(19)

Let

$${\mathit{k}}_{\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}={\int }_{A_{E0}}{\mathit{B}}^{\mathrm{T}}\left(x,y\right) {\mathit{D}}_0 \mathit{T} \mathrm{d}{A}_{E0}.$$

(20)

Then ${\mathit{k}}_{nonlocal}$ can be interpreted as the element coupling stiffness contribution (based on the center-point approximation) used to assemble the nonlocal global stiffness matrix. For the entire plate structure, the nonlocal global stiffness matrix ${\mathit{K}}_{\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}$ is assembled from all element coupling stiffness contributions and linearly combined with the local global stiffness matrix ${\mathit{K}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}$ according to the Eringen two-phase model

$$\mathit{K}=\left(1-{\xi }_2\right) {\mathit{K}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}+{\xi }_2 {\mathit{K}}_{\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l}}.$$

(21)

The model under consideration naturally degenerates into the classical local Kirchhoff plate finite element when the dimensionless nonlocal parameter, denoted by $\kappa \to 0$, or when ${\xi }_2\to 0$. while still simultaneously reflecting the long-range coupling effects between elements under finite nonlocal parameters.

It should be noted that the above derivation is based on the element center point approximation and equal area assumption, allowing the nonlocal integral constitutive to be embedded in the standard Kirchhoff plate finite element framework in the form of a matrix.

3.3. Nonlocal Kernel Functions and the Selection of Their Influence Domains

In the nonlocal integral constitutive relation (Equation 1) and the aforementioned finite element discretization format (Equation 11), the following kernel function $\alpha \left(|x^{\prime }-x|,\kappa \right)$ plays a decisive role. Physically, it describes the nonlocal weighted influence of the strain field at the source point $x^{\prime }$ on the stress state at the reference point $x$. To ensure the physical completeness of the nonlocal theory, the kernel function must satisfy three basic properties: The kernel function reaches a maximum at $x^{\prime }=x$ and decays monotonically as the Euclidean distance $s=|x^{\prime }-x|$ between the two points increases. This reflects the physical fact that long-range interactions between microscopic particles weaken with increasing distance; when the dimensionless nonlocal parameter $\kappa \to 0$, the kernel function should degenerate into the Dirac $\delta$ function, at which point the nonlocal elasticity theory reverts to the classical local elasticity theory. Therefore, classical elasticity theory can also be viewed as a special case of the nonlocal theory when long-range forces are ignored; to ensure the completeness of the constitutive relation, the integral of the kernel function over the entire domain should be 1:

$${\int }_V\alpha \left(\right|x^{\prime }-x|,\kappa )\mathrm{d}V\left(x^{\prime }\right)=1.$$

(22)

While any function that satisfies the aforementioned three properties can theoretically serve as a kernel function, it is important to note that different kernel function forms result in different constitutive responses. In this paper, the two most commonly used two-dimensional isotropic kernel functions are examined: the exponential and Gaussian kernels. These functions are selected for analysis of the two-dimensional Kirchhoff plate problem. The mathematical expressions pertaining to these phenomena are as

$$\alpha \left(s,l_a\right)={\displaystyle \frac{1}{2\pi l_a^2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{s}{l_a}}\right),$$

(23)

$$\alpha \left(s,l_a\right)={\displaystyle \frac{1}{4\pi l_a^2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{s^2}{4l_a^2}}\right),$$

(24)

where $l_a=e_0{a}_0$ is the nonlocal parameter.

In numerical implementation, considering the decay characteristic of the kernel function, when the distance $s$ exceeds a specific range, its function value rapidly approaches zero. To reduce the computational cost, we truncate the kernel by introducing a finite influence domain. After truncation, to satisfy equation (22), it must be normalized again. This paper sets the two-dimensional nonlocal influence domain as a circular region centered at the reference point with a radius $R=3l_a$. Nonlocal interactions outside this region are considered negligible. When taking $l_a=1\mathrm{n}\mathrm{m}$, the three-dimensional plots of the above two kernel functions are shown in Figure 1 and Figure 2, respectively.

4. Numerical Results

In the preceding sections, a nonlocal Kirchhoff plate finite element model was established based on inter-belt analysis. This model was used to derive the stiffness matrix and governing equations. To verify the correctness, convergence, and effectiveness of this algorithm in handling nonlocal effects, a series of numerical examples are presented in this section.

In the numerical examples of this section, unless otherwise stated, the two-dimensional exponential kernel function is selected. The geometric and material parameters of the plate are set as follows: elastic modulus $E=1.06\times {10}^{12}\mathrm{Pa}$, Poisson's ratio $\nu =0.25$, mass density $\rho =2250{\mathrm{kg}/\mathrm{m}}^3$. To consider the small-scale effect, the nonlocal parameter is $l_a=e_0{a}_0$, and the two-phase mixing parameter is ${\xi }_2$. The frequency parameter is defined as ${\omega }_0=\omega a^2\sqrt{\rho h/D_0}$, where $a=10\mathrm{n}\mathrm{m}$ is the width of the square plate, $h=0.34\mathrm{n}\mathrm{m}$ is the plate thickness, and $D_0$ is the bending stiffness of the plate.

4.1. Algorithm Convergence and Accuracy Verification

To verify the accuracy of the nonlocal finite element algorithm based on inter-belt analysis proposed in this paper, the following section first examines the convergence behavior of the frequency parameters under different mesh densities, then compares the calculation results with the existing literature, and determines the optimal two-phase mixing parameter for this kernel function.

As demonstrated in Figure 3, the frequency parameter of the Kirchhoff plates converges with increasing element number under simply supported boundary conditions on all four sides (SSSS). This is achieved using the exponential kernel function and varying nonlocal parameter values. Figure 3(a) illustrates the trend of the first-order frequency parameter. The variables $m$ and $n$ denote the number of half-waves of the plate along the x and y axes, respectively. Their magnitude determines the magnitude of the mode index. The three curves in the figure correspond to the local case, i.e. when the nonlocal parameter $l_a=0\mathrm{nm}$, and cases where the nonlocal parameter $l_a=1\mathrm{nm}$ and $2\mathrm{nm}$, respectively. The mesh division increases sequentially from 10×10 to 60×60.

Results indicate that, under local conditions, as the number of elements increases, the calculated first-order frequency parameter exhibits a gradual increase and stabilizes after 30×30 elements, indicating good convergence. Following the introduction of the nonlocal effect, a modification in the convergence behavior of the frequency parameter becomes apparent. For the nonlocal parameter equal to 1 nm, the frequency parameter decreases gradually as the element number increases. This decrease stabilizes once the element number reaches 40×40. Furthermore, when the nonlocal parameter is increased to 2 nm, the frequency parameter exhibits a similar trend of decreasing with the increase of element number. However, the convergence speed is slightly faster, stabilizing after the element number reaches 30×30. This indicates that incorporating nonlocal effects not only modifies the dynamic characteristics of the structure but also significantly affects the convergence of numerical calculations. However, the algorithm presented in this paper can attain stable solutions under suitable mesh densities.

Figure 3(b) provides a detailed illustration of the convergence of the fourth-order frequency parameter when other conditions remain constant. The findings indicate that, under the examined nonlocal parameters, the fourth-order frequency parameters can also converge to stable values, once the element number reaches 30×30. It is evident from the convergence analysis that, unless stated otherwise, a 30×30 mesh division is to be employed in subsequent numerical validations. This ensures calculation accuracy and efficiency.

The verification of algorithm convergence serves as a foundation for the subsequent assessment of algorithmic accuracy. As illustrated in Figure 4, the trend of the first-order frequency parameter with the nonlocal parameter is demonstrated under different two-phase mixing parameters, namely, ${\xi }_2$. This is then compared with the finite element results obtained by Shahidi et al. (2013) [22]. The selection of ${\xi }_2$ follows literature recommendations and is used to calibrate the consistency between integral two-phase models and commonly used differential models.

As demonstrated in Figure 4(a), with the alteration of the two-phase mixing parameter, there is a substantial discrepancy in the frequency parameter curves. This figure shows that when the value of ${\xi }_2$ falls between 0.55 and 0.56, the curve calculated in this paper exhibits the highest degree of agreement with the data presented in the literature. To achieve a more precise quantification of this agreement, Figure 4(b) presents a further analysis of the average relative error between the calculated results and the literature’s different two-phase mixing parameters. The findings suggest that the average relative error attains its minimum when the parameter is approximately 0.56, at which point the average relative error is of the order of ${10}^{-3}$. This result provides a significant reference point for the subsequent determination of the optimal two-phase mixing parameter for this kernel function.

It is noteworthy that, despite the theoretical value range of the two-phase mixing parameter ${\xi }_2$ being $\left[\mathrm{0,1}\right]$, as${\xi }_2\to 1$, the model degenerates to the pure integral nonlocal formulation. This theoretical framework has been demonstrated to give rise to ill-posed problems within bounded domains. Consequently, to ensure physical consistency and numerical stability of the results, it is recommended in related studies to set ${\xi }_2$ around 0.5 [46,47]. Using the exponential kernel, the parameter is calibrated against the benchmark by minimizing the average relative error and is approximately selected as ${\xi }_2=0.56$, consistent with the literature-recommended range around 0.5.

As illustrated in Figure 5, the calculation results are presented when using the Gaussian kernel function. In a manner analogous to the exponential kernel function, the frequency parameters under the Gaussian kernel function are also contingent on the two-phase mixing parameter ${\xi }_2$. Error analysis shows that, under this kernel function, when approximately 0.52 is attained by ${\xi }_2$, the relative error between the present frequency parameters and those reported in the literature is minimized. Despite the fact that the optimal values of ${\xi }_2$ corresponding to different kernel functions exhibit slight discrepancies (0.56 for the exponential kernel, 0.52 for the Gaussian kernel), both are situated in proximity to 0.5, thereby substantiating the rationality of this parameter value. Despite the superior smoothness of the Gaussian kernel function in physical contexts, the exponential kernel function occupies a more central position in nonlocal theory due to its distinctive mathematical properties. This is because the exponential kernel function can accurately fit the dispersion curves of atomic lattice dynamics. Moreover, the Green’s function of a linear differential operator provides mathematical equivalence between the nonlocal integral constitutive relation and high-order differential equations [48]. This characteristic enables the transformation of originally complex integral equations into differential forms, which are more amenable to solution. This, in turn, greatly facilitates theoretical analysis and the acquisition of analytical solutions. In contrast, the Gaussian kernel function is unable to provide this mathematical bridge for integral-differential transformation and typically relies exclusively on numerical solutions. Therefore, to facilitate comparison with existing differential-form nonlocal studies, the subsequent parametric analysis is mainly based on the exponential kernel function, with ${\xi }_2=0.56$ to obtain relatively accurate solutions.

As demonstrated in the preceding convergence analysis and the subsequent determination of the optimal two-phase mixing parameter, the tables (1 to 4) present a comparison of the frequency parameters calculated in this study with existing literature results under the exponential kernel function (Table 1 and Table 2) and the Gaussian kernel function (Table 3 and Table 4), respectively. The present study compares the relative errors of the results obtained in this paper with those reported in two other papers: Shahidi et al. (2013) [22] and Pradhan and Phadikar (2009) [10].

From Table 1, Table 2, Table 3 and Table 4, it can be seen that based on the two-phase mixing parameter determined previously (${\xi }_2=0.56$ for exponential kernel, ${\xi }_2=0.52$ for Gaussian kernel), the frequency parameters calculated using the inter-belt analysis method in this paper agree well with the results in the literature, with relative deviations in most cases being below 2%. In the local case where the nonlocal parameter is 0, the relative error is only of the order of ${10}^{-3}$, thereby verifying the accuracy of the algorithm's degeneration. Further observations of the error distribution under different nonlocal parameters revealed that both kernel functions showed remarkably high agreement when the nonlocal parameter $l_a=2$. This was evidenced by relative errors falling below 0.1%, thereby demonstrating the model's high consistency under this particular parameter. In the context of larger nonlocal parameters ($l_a=\mathrm{3,4}$), the relative error of the exponential kernel function results (approximately 0.4% - 0.8%) is marginally lower than that of the Gaussian kernel function (approximately 0.8% - 0.9%). This finding suggests that the exponential kernel function is numerically more stable in this range. Overall, for the two kernel functions considered and the tested range of nonlocal parameters, the relative errors remain within a reasonable range. This confirms the accuracy and reliability of the proposed integral-form nonlocal finite element algorithm for the two kernel functions considered in this study and the nonlocal parameters examined.

4.2. Parameter Study

4.2.1. Influence of Nonlocal Parameters and Two-Phase Mixing Parameter

As demonstrated in Figure 6, the variation of the percentage reduction of the first five frequency parameters with the nonlocal parameter $l_a$ is observed under the four-sided simply supported (SSSS) boundary conditions. These results show that as the nonlocal parameter increases, the percentage reduction in frequency for each mode exhibits a marked upward trend. This finding suggests that enhancing the nonlocal effect will reduce the plate's overall stiffness, thereby decreasing its natural frequency. Furthermore, the frequency reduction amplitude of higher-order modes is significantly larger than that of lower-order modes, indicating that high-frequency vibration is more sensitive to nonlocal effects, which is consistent with the general principle of nonlocal elasticity that short-wavelength deformation is more affected by micro-scale effects.

Figure 7 further illustrates the effect of the two-phase mixing parameter ${\xi }_2$ on the percentage reduction in the frequency parameters of the initial five vibration modes when the nonlocal parameter, $l_a=1\mathrm{n}\mathrm{m}$, is constant. The results demonstrate that as the parameter ${\xi }_2$ increases from 0.1 to 0.9, the frequency-reduction percentage increases linearly. As the nonlocal integral term is represented by the variable ${\xi }_2$, an increase in its value directly strengthens the nonlocal softening effect, thereby leading to a further reduction in structural stiffness and vibration frequency.

4.2.2. Influence of Aspect Ratio

As illustrated in Figure 8, the variation curve of the first-order frequency parameter of the plate with the aspect ratio is depicted for different sizes of nanoplates under a fixed nonlocal parameter $l_a=1\mathrm{n}\mathrm{m}$. In this analysis, the plate thickness and nonlocal parameter are kept constant. It has been observed that as the aspect ratio increases, the first-order natural frequency decreases monotonically. When the aspect ratio is small, the frequency declines precipitously (as illustrated by the steeper curve in the figure); conversely, when the aspect ratio exceeds a certain threshold, the frequency change curve tends to flatten. This phenomenon shows that, at the nanoscale, alterations in the plate's geometric configuration exert a substantial influence on its dynamic characteristics. However, as the plate becomes more slender, the marginal impact of this geometric effect diminishes concomitantly. It is notable that when the aspect ratio exceeds 5, the nonlocal plate's influence on the frequency becomes less significant.

4.2.3. Influence of Thickness Ratio

As illustrated in Figure 9, the thickness ratio significantly influences the first-order frequency parameter of nanoplates of varying dimensions, with the nonlocal parameter fixed at 1 nm. The thickness ratio range for the numerical example is set to $\left[0.034,0.1\right]$. The lower limit corresponds to the benchmark thickness parameter for single-layer graphene, with the aim of verifying the model's effectiveness for single-atomic-layer nanomaterials. The upper limit corresponds to the upper limit of the thin plate applicable range, which is usually defined for Kirchhoff thin plates ($h/a\le 1/10$). This ensures the rationality of the assumption of neglecting transverse shear deformation. The results show that, as the thickness ratio increases, the first-order natural frequency of the plate increases approximately linearly. This is due to the fact that an increase in thickness directly enhances the bending stiffness of the plate, thereby increasing the natural frequency of the structure. It is noteworthy that in the presence of nonlocal effects, the influence of thickness ratio variation on frequency still adheres to the fundamental principles of classical mechanics. That is to say, stiffness enhancement results in an increase in frequency.

5. Discussion

This study addresses the computational difficulty of Eringen’s integral-form nonlocal elasticity for plate problems, where the governing equations are integro-differential and require cumbersome global integration. A finite element formulation based on symplectic inter-belt discretization is proposed for the free-vibration analysis of two-dimensional nonlocal Kirchhoff plates. By discretizing long-range interactions as element-to-element coupling within a prescribed influence domain, the nonlocal integral operator is embedded into a standard finite element framework, allowing the resulting governing equations to be cast into a conventional linear algebraic eigenvalue problem. The formulation is applicable to different rational kernel functions, nonlocal parameters, and two-phase mixture parameters, and thus provides a numerical approach for integral nonlocal vibration analysis of nanoplates.

The convergence and accuracy of the proposed method are verified through mesh refinement studies and comparisons with published results. The numerical results consistently indicate that nonlocal interactions lead to stiffness softening, manifested by a reduction in natural frequencies compared with the local model. Moreover, higher-order modes exhibit stronger sensitivity to the nonlocal parameter, which aligns with the general observation in nonlocal mechanics that short-wavelength responses are more affected by small-scale effects. The parametric investigations further clarify how the nonlocal parameter, mixture parameter, kernel-function type, and geometric ratios jointly influence the vibration characteristics, providing quantitative guidance for selecting modeling parameters in practical analyses.

Beyond the methodological contribution, the proposed framework is relevant to nanoengineering applications where resonant frequencies must be predicted with high fidelity. In MEMS/NEMS devices such as nanomechanical resonators, mass sensors, and RF filters, small deviations in stiffness can cause noticeable frequency shifts, especially at high frequencies. Classical local continuum models may overestimate the effective stiffness of nanostructures by neglecting long-range interactions. In contrast, the present integral nonlocal formulation can provide frequency predictions that reflect size-dependent softening, thereby supporting frequency calibration, device sizing, and material selection for nanoscale components.

Several limitations of the present work should be acknowledged. The formulation is developed under linear elasticity and is demonstrated primarily for regular rectangular geometries. Geometric nonlinearity associated with large deflection, as well as multi-physics couplings and more complex boundary conditions, are not considered here and may be important for certain operating regimes and device designs. Future work will extend the inter-belt framework to nonlinear vibration and multi-physics problems and will explore integration with advanced discretization techniques (e.g., isogeometric analysis) to better handle complex geometries and material heterogeneity.

6. Conclusions

A symplectic inter-belt finite element formulation is developed for free-vibration analysis of Kirchhoff plates governed by Eringen’s integral nonlocal elasticity. By discretizing long-range interactions as element-to-element couplings within a finite influence domain, the integral operator is incorporated into a standard finite element framework, yielding a generalized eigenvalue problem for modal analysis.

Mesh-refinement studies show stable convergence of the frequency parameters for the investigated cases. With calibrated two-phase mixing parameters, benchmark comparisons indicate that the present results agree well with published solutions: relative deviations are below $2\%$ in most cases, and the local limit $(l_a=0)$yields errors on the order of ${10}^{-3}$. The calibration based on average relative error suggests ${\xi }_2\approx 0.56$ for the exponential kernel and ${\xi }_2\approx 0.52$ for the Gaussian kernel, both close to the literature-recommended range around 0.5.

Parametric results confirm that stiffness softening and frequency reduction occur as the nonlocal length scale increases, with higher-order modes showing greater sensitivity to nonlocal effects. The frequency-reduction percentage increases approximately linearly with ${\xi }_2$ (from 0.1 to 0.9). At the same time, geometric studies show a monotonic decrease of the first-order frequency with increasing aspect ratio (with reduced sensitivity beyond an aspect ratio of about 5) and an approximately linear increase with thickness ratio over $h/a\in [0.034,0.1]$.