**1. Introduction**

As typical structures, thick plates with holes are widely used under dynamic loads in aerospace, ocean engineering, civil engineering and mechanical engineering [1–5]. Due to the stress concentration near the load-bearing opening in the finite thickness structure, there will be an intense three-dimensional effect zone. The three-dimensional effect zone is closely related to the relative thickness of the structure, which largely controls the fracture, fatigue and other mechanical properties of the structure [6,7]. There is a large error in calculating the dynamic problem of the plate with actual thickness based on the classical thin plate theory. With the development of modern science and technology, engineering structure designs tend to be light, and the way to achieve this is to use advanced materials and improve the structure design theory. However, because the three-dimensional problem has not been well solved, it will encounter great mathematical difficulties in solving the three-dimensional problem.

The classical plate theory (CPT) is the basic theory in the hierarchy of plate theories [6]. Since it was formulated systematically in the 19th century, CPT has been widely applied for buckling, bending and vibration analyses of plates. However, it is worth noting that CPT does not consider the effects of shear deformation or rotary inertia [8–10]. Hence, in the analysis of thick plates and also in the case of thin plates vibrating at higher frequencies, the use of CPT would result in considerable errors [11–14]. These limitations of CPT led to the development of first-order and higher-order shear deformation plate theories.

Reissner [11,12] proposed a thick plate theory, including the effects of shear deformation. Reissner's theory involves three coupled governing differential equations in terms of

**Citation:** Zhou, C.; Wang, M.; Han, X.; Xue, H.; Ni, J.; Zhou, W. A Novel Exact Plate Theory for Bending Vibrations Based on the Partial Differential Operator Theory. *Mathematics* **2021**, *9*, 1920. https:// doi.org/10.3390/math9161920

Academic Editor: Francisco Beltran-Carbajal

Received: 24 May 2021 Accepted: 2 August 2021 Published: 12 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

19

three unknowns. In fact, Reissner's theory is based on a stress approach. Mindlin [13] considered both the effects of shear deformation and rotary inertia and proposed a first-order displacement-based theory, which involves three coupled governing differential equations. Mindlin's theory was derived by using the frequency domain method, so its disadvantage is that it cannot predict the constant transverse shear strains and stresses across the plate thickness. Therefore, the theory requires a shear correction factor to match the strain energy calculated by using constant transverse shear strains and stresses. The transverse shear strains and stresses vary parabolically across the plate thickness [14].

The Levinson [15] and Reddy [16,17] plate theories are higher-order shear deformation plate theories based on displacement. Both the theories are based on the same displacement field. However, the governing equations of the Levinson theory are derived from the plate gross equilibrium equations, whereas Reddy used the principle of virtual displacements to derive the governing differential equations. Neither theory requires a shear correction factor. The plate theory proposed by Kant [18] is based on a higher-order displacement model, which causes a secondary change in transverse shear strain and a linear change in transverse normal strain across the plate thickness. This work involves the flexure of thick rectangular isotropic plates. The formulation of the theory involves six variables. Since the higher-order plate theories cannot take the continuity conditions of displacement and shear stresses into consideration, the accuracy of analytical solutions cannot meet the requirements in engineering applications [19].

Nedri [20] presented a novel refined hyperbolic shear deformation theory based on the assumption that the transverse displacements consist of shear and bending components where the bending components do not contribute to shear forces, and likewise, the shear components do not contribute to bending moments. The irrelevancies of the two components make calculations simple, yet it may lead to an erroneous result when the displacements vary sharply across the thickness. Shimpi [21] developed a refined plate theory (RPT) with only two unknown functions available in the paper. RPT produces two fourth-order governing differential equations, which are uncoupled for static problems and are only inertially coupled for dynamic problems. Subsequently, Shimpi [22] presented two new first-order shear deformation plate theories with only two unknown functions to improve RPT. Using the method proposed by Timoshenko and Ashwell, Nicassio [4] presented a novel forecast model to map the surface profiles of bistable laminates and developed an analytical model to provide an interpretation of the bistable shapes in terms of principal and anticlastic curvatures. Wu [5] gave a revised method to increase the stiffness and natural frequency of bistable composite shells, which can be suitable for spatial, lightweight structural components.

Among various refined plate theories, Carrera [23] presented a unified expression, which provides a program to obtain refined structural models of beams, plates and shells that explain variable kinematics descriptions. Based on Carrera's unified formulation (CUF), these structural models are obtained using the N-order Taylor expansion to expand the unknown displacement variables. Tornabene et al. [24,25] derived a general formulation of 2D higher-order equivalent plate theory. The theoretical framework covers the static and dynamic analysis of shell structures by using a general displacement field based on CUF. In the CUF system, a linear case can describe a classical model, while a higherorder case can describe a three-dimensional structure. Kolahchi [3] investigated bending, buckling and buckling of embedded nano-sandwich plates based on refined zigzag theory (RZT), sinusoidal shear deformation theory (SSDT), first-order shear deformation theory (FSDT) and classical plate theory (CPT), and a differential cubature (DC) method is applied for obtaining the static response, the natural frequencies and the buckling loads of nano-sandwich plates. The numerical investigation shows that RZT is highly accurate in predicting the deflection, frequency and buckling load of nano-sandwich plates without requiring any shear correction factors.

In brief, although the above CPT theories have been optimized and updated, they are theoretically based on the geometric method, and the models are rough since engineering

assumptions are still used during the derivation, which results in many limitations in the application of thick wall structures, especially in the case of plates vibrating at higher frequencies. In this paper, we propose a novel theory of exact elastic dynamics for bending plates not based on the geometric view but the algebraic view. During the derivation, we apply the general formal solution proposed by Boussinesq–Galerkin and the operator theory of partial differential equations. The exact elasto-dynamics equations for bending plates are obtained by using appropriate gauge conditions, and the exact dynamic theory of thick plates is compared with other plate theories. Since the derivation of the dynamic equation is carried out without any prior assumption, the proposed dynamic equation of plates is more exact and can be applied in a wider frequency range and greater thickness. The exact thick plate theory in this paper makes up for the shortcomings of the classical thin plate theory and other thick plate theories. It can be used not only for structures with large thickness span ratio but also for vibration mechanics problems with great influence of shear deformation and moment of inertia, such as spacecraft attitude dynamics and control, structural motion stability, dynamic stress concentration for thick plates with holes, and calculations for submarine anechoic tile structure design.

#### **2. Derivation of the Exact Dynamic Theory for the Bending Plate**

First, the derivation process of the exact plate theory for the bending plate is introduced. According to the three-dimensional elasto-dynamics theory, the governing equation of the spatial displacement field is the Navier equation:

$$
\mu \nabla^2 \mathbf{u} + (\lambda + \mu) \nabla(\nabla \cdot \mathbf{u}) = \rho \frac{\partial^2 \mathbf{u}}{\partial t^2},
\tag{1}
$$

where *μ*, *λ* are the Lame constants, ∇ = **i***∂*/*∂x* + **j***∂*/*∂y* + **k***∂*/*∂z*, *ρ* is the density.

From Equation (1), based on Boussinesq–Galerkin solution (B-G solution), the solution given as:

$$\mathbf{u} = 2(1 - \nu) \left( \nabla^2 - \frac{1}{c\_1^2} \frac{\partial^2}{\partial t^2} \right) \mathbf{G} - \nabla(\nabla \cdot \mathbf{G}),\tag{2}$$

where *c*1, *c*<sup>2</sup> are longitudinal wave velocity and transverse wave velocity, *ν* is the Poisson ratio, and **G** = (*G*1, *G*2, *G*3) is the Somigliana vector potential function, which satisfies the following relation as:

$$\prod \left(\nabla^2 - T\_{\vec{\gamma}}^2\right) \mathbf{G} = 0. \tag{3}$$

where *Tj* are time differential operators, *T*<sup>2</sup> *<sup>j</sup>* <sup>=</sup> <sup>1</sup> *c*2 *j ∂*2 *<sup>∂</sup>t*<sup>2</sup> ,(*<sup>j</sup>* = 1, 2),

Using the Taylor series expansion of the exponential operator function, the displacement at any point in the plate can be written as:

$$
\mu\_x(x, y, z) = \exp\left(z \frac{\partial}{\partial z}\right) \mu\_x(x, y, 0), \tag{4}
$$

$$
\mu\_y(x, y, z) = \exp\left(z \frac{\partial}{\partial z}\right) \mu\_y(x, y, 0),
\tag{5}
$$

$$u\_{\overline{z}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \exp\left(z \frac{\partial}{\partial z}\right) u\_{\overline{z}}(\mathbf{x}, \mathbf{y}, \mathbf{0}). \tag{6}$$

The fluctuation of plate bending is a case of antisymmetric motion, and Equation (4) can be written as:

$$
\mu\_x(x, y, z) = \sin \hbar \left( z \frac{\partial}{\partial z} \right) \mu\_x(x, y, 0),
\tag{7}
$$

$$
\mu\_{\mathcal{Y}}(x, y, z) = \sin \mathbf{h} \left( z \frac{\partial}{\partial z} \right) \mu\_{\mathcal{Y}}(x, y, 0), \tag{8}
$$

$$u\_z(x,y,z) = \cosh\left(z\frac{\partial}{\partial z}\right)u\_z(x,y,0),\tag{9}$$

where sinh(·) is hyperbolic sine function, and cosh(·) is a hyperbolic cosine function.

The B-G solution can be written as:

$$\mathbf{G}\_{\dot{\jmath}}(\mathbf{x}, y, \mathbf{z}) = \sin \mathbf{h} \left( z \frac{\partial}{\partial z} \right) \mathbf{G}\_{\dot{\jmath}}(\mathbf{x}, y, 0) = \cosh \left( z \frac{\partial}{\partial z} \right) \sum\_{i=1}^{2} \mathbf{G}\_{\dot{\jmath}}^{i}(\mathbf{x}, y, 0), \tag{10}$$

$$\mathbf{G}\_3(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \cosh\left(z \frac{\partial}{\partial z}\right) \mathbf{G}\_3(\mathbf{x}, \mathbf{y}, \mathbf{0}) = \cosh\left(z \frac{\partial}{\partial z}\right) \sum\_{i=1}^2 \mathbf{G}\_3^i(\mathbf{x}, \mathbf{y}, \mathbf{0}),\tag{11}$$

The Somigliana vector potential function **G** can be decompose into two vector potential function as **<sup>G</sup>** <sup>=</sup> **<sup>G</sup>**<sup>1</sup> <sup>+</sup> **<sup>G</sup>**2, and ∇2 *<sup>j</sup>* <sup>+</sup> *<sup>∂</sup>*<sup>2</sup> *∂z*<sup>2</sup> **<sup>G</sup>***<sup>j</sup>* <sup>=</sup> 0, here <sup>∇</sup><sup>2</sup> *<sup>j</sup>* <sup>=</sup> <sup>∇</sup><sup>2</sup> <sup>−</sup> *<sup>T</sup>*<sup>2</sup> *<sup>j</sup>* ,(*j* = 1, 2) is the Lorentz operator. The trigonometric function operator can be written as:

$$\frac{\sin\left(z\nabla\_{\dot{j}}\right)}{\nabla\_{\dot{j}}} = \sum\_{n=0}^{\infty} (-1)^{n} \frac{1}{(2n+1)!} z^{2n+1} \nabla\_{\dot{j}}^{2n} \, , \tag{12}$$

$$\cos\left(z\nabla\_{\dot{\slash}}\right) = \sum\_{n=0}^{\infty} (-1)^n \frac{1}{(2n)!} z^{2n} \nabla\_{\dot{\slash}}^{2n},\tag{13}$$

here *j* = 1, 2.

We can also obtain the following relation as:

$$\nabla\_0 \mathbf{G} = \frac{\partial G\_1}{\partial \mathbf{x}} + \frac{\partial G\_2}{\partial y} + \frac{\partial G\_3}{\partial z} = \sum\_{j=1}^2 \frac{\sin \left( z \nabla\_j \right)}{\nabla\_j} \left( \frac{\partial \mathbf{g}\_1^j}{\partial \mathbf{x}} + \frac{\partial \mathbf{g}\_2^j}{\partial y} - \nabla\_j^2 \mathbf{g}\_3^j \right). \tag{14}$$

For the sake of avoiding the non-uniqueness of unknown functions, two gauge conditions are adopted as follows:

$$\frac{\partial \mathcal{g}\_1^j}{\partial x} + \frac{\partial \mathcal{g}\_2^j}{\partial y} = 0,\\
(j = 1, 2). \tag{15}$$

Equation (8) can be written as:

$$\begin{array}{rcl} \nabla\_0(\nabla\_0 \cdot \mathbf{G}) &=& -\sum\_{j=1}^2 \left[ \frac{\sin(z \nabla\_j)}{\nabla\_j} \nabla\_j^2 \frac{\partial}{\partial x} g\_3^j \right] \mathbf{i} - \sum\_{j=1}^2 \left[ \frac{\sin(z \nabla\_j)}{\nabla\_j} \nabla\_j^2 \frac{\partial}{\partial y} g\_3^j \right] \mathbf{j} \\ &- \sum\_{j=1}^2 \left[ \cos(z \nabla\_j) \nabla\_j^2 g\_3^j \right] \mathbf{k} \end{array} \tag{16}$$

The displacement in the plate can be expressed as:

$$\mathbf{u} = 2(1 - \nu) \left( \nabla\_1^2 + \frac{\partial^2}{\partial z^2} \right) \mathbf{G} - \nabla\_0 (\nabla\_0 \cdot \mathbf{G}).\tag{17}$$

Its component-wise expressions can be written as:

$$\ln \mu\_{\mathbf{x}} = \frac{\sin(z \nabla\_{2})}{\nabla\_{2}} T\_{2}^{2} g\_{1}^{2} + \sum\_{j=1}^{2} \frac{\sin(z \nabla\_{j})}{\nabla\_{j}} \nabla\_{j}^{2} \frac{\partial}{\partial \mathbf{x}} g\_{3 \prime}^{j} \tag{18}$$

$$u\_y = \frac{\sin(z\nabla\_2)}{\nabla\_2} T\_2^2 g\_2^2 + \sum\_{j=1}^2 \frac{\sin(z\nabla\_j)}{\nabla\_j} \nabla\_j^2 \frac{\partial}{\partial y} g\_{3'}^j \tag{19}$$

$$u\_z = \cos(z\nabla\_2) T\_2^2 g\_3^2 + \sum\_{j=1}^2 \cos(z\nabla\_j) \nabla\_j^2 g\_3^j. \tag{20}$$

Considering the neutral surface displacement and the normal angle, the generalized displacement in the plate can be expressed as:

$$\psi\_{\mathbf{x}} = \left. -\frac{\partial u\_{\mathbf{x}}}{\partial z} \right|\_{z=0} = -T\_2^2 g\_1^2 - \sum\_{j=1}^2 \nabla\_j^2 \frac{\partial^2 g\_3^j}{\partial \mathbf{x}},\tag{21}$$

$$\psi\_y = -\frac{\partial u\_y}{\partial z}\Big|\_{z=0} = -T\_2^2 \text{g}\_2^2 - \sum\_{j=1}^2 \nabla\_j^2 \frac{\partial^2 \text{g}\_3^j}{\partial y},\tag{22}$$

$$w = \left. u\_z \right|\_{z=0} = \nabla^2 \mathfrak{g}\_3^2 + \nabla\_1^2 \mathfrak{g}\_3^1. \tag{23}$$

The rotational normal angle to the neutral surface can be expressed as:

$$
\psi\_{\mathbf{x}} = \frac{\partial F}{\partial \mathbf{x}} + \frac{\partial f}{\partial y}, \quad \psi\_y = \frac{\partial F}{\partial y} - \frac{\partial f}{\partial \mathbf{x}}.\tag{24}
$$

The functions *g*<sup>2</sup> <sup>1</sup>, *<sup>g</sup>*<sup>2</sup> <sup>2</sup>, *<sup>g</sup>*<sup>2</sup> <sup>3</sup> can be expressed by the neutral surface displacement and normal angle as:

$$g\_1^2 = -T\_2^{-2} \left(\frac{\partial f}{\partial y} + \frac{\partial E}{\partial \mathbf{x}}\right),\tag{25}$$

$$
\delta g\_2^2 = T\_2^{-2} \left( \frac{\partial f}{\partial x} + \frac{\partial E}{\partial y} \right) \tag{26}
$$

$$g\_3^1 = -T\_2^{-2} \nabla\_1^{-2} \left(\nabla\_1^2 w + \nabla^2 F\right),\tag{27}$$

$$g\_3^2 = T\_2^{-2}(F + w - E). \tag{28}$$

where <sup>∇</sup>2*<sup>E</sup>* <sup>=</sup> 0,*<sup>F</sup>* <sup>=</sup> −∇<sup>2</sup> 1*g*1 <sup>3</sup> − ∇<sup>2</sup> 2*g*2 <sup>3</sup> + *E*.

In this way, the displacement can be derived as:

$$\begin{split} u\_{X} &= \sum\_{j=1}^{2} (-1)^{j-1} \frac{\sin(z\nabla\_{j})}{\nabla\_{j}} \frac{\partial w}{\partial x} - \frac{\sin(z\nabla\_{2})}{\nabla\_{2}} \left( \frac{\partial F}{\partial x} + \frac{\partial w}{\partial x} \right) \\ &- T\_{2}^{-2} \sum\_{j=1}^{2} \left( -1 \right)^{j-1} \frac{\sin(z\nabla\_{j})}{\nabla\_{j}} \nabla^{2} \left( \frac{\partial F}{\partial x} + \frac{\partial w}{\partial x} \right) \end{split} \tag{29}$$

$$\begin{split} u\_{\mathcal{Y}} &= \sum\_{j=1}^{2} (-1)^{j-1} \frac{\sin(z \nabla\_{j})}{\nabla\_{j}} \frac{\partial w}{\partial \mathcal{Y}} - \frac{\sin(z \nabla\_{2})}{\nabla\_{2}} \left( \frac{\partial F}{\partial \mathcal{Y}} - \frac{\partial f}{\partial x} \right) \\ &- T\_{2}^{-2} \sum\_{j=1}^{2} (-1)^{j-1} \frac{\sin(z \nabla\_{j})}{\nabla\_{j}} \nabla^{2} \left( \frac{\partial F}{\partial y} + \frac{\partial w}{\partial y} \right) \end{split} \tag{30}$$

$$u\_z = \cos(z\nabla\_1)w - T\_2^{-2}\sum\_{j=1}^2 (-1)^{j-1} \cos(z\nabla\_2)\nabla^2(F+w). \tag{31}$$

In line with Hooke's law, the stress components can be expressed as:

$$\begin{split} \tau\_{zx} &= 2\mu \cos(z\nabla\_1) \frac{\partial w}{\partial x} - \mu \cos(z\nabla\_2) \left( \frac{\partial F}{\partial x} + \frac{\partial w}{\partial x} + \frac{\partial f}{\partial y} \right) \\ &- 2\mu T\_2^{-2} \sum\_{j=1}^2 (-1)^{j-1} \cos(z\nabla\_j) \nabla^2 \left( \frac{\partial F}{\partial x} + \frac{\partial w}{\partial x} \right) \end{split} \tag{32}$$

$$\begin{split} \pi\_{\overline{z}\underline{y}} &= 2\mu \cos(z\nabla\_{1}) \frac{\partial \underline{w}}{\partial \underline{y}} - \mu \cos(z\nabla\_{2}) \left( \frac{\partial F}{\partial \underline{y}} + \frac{\partial \underline{w}}{\partial \underline{y}} - \frac{\partial f}{\partial \underline{x}} \right) \\ &- 2\mu T\_{2}^{-2} \sum\_{j=1}^{2} (-1)^{j-1} \cos(z\nabla\_{j}) \nabla^{2} \left( \frac{\partial F}{\partial \underline{y}} + \frac{\partial \underline{w}}{\partial \underline{y}} \right) \end{split} \tag{33}$$

$$\begin{cases} \sigma\_z = 2\mu T\_2^{-2} \sum\_{j=1}^2 (-1)^{j-1} \frac{\sin(z\nabla\_j)}{\nabla\_j} \nabla^2 \nabla^2 (F + w) \\ \quad - (\lambda + 2\mu) T\_1^{-2} T\_2^{-2} \frac{\sin(z\nabla\_1)}{\nabla\_1} \nabla^2 (F + w) \\ \quad + (\lambda + 2\mu) T\_1^2 \frac{\sin(z\nabla\_1)}{\nabla\_1} w \\ \quad - 2\mu \frac{\sin(z\nabla\_1)}{\nabla\_1} \nabla^2 w + 2\mu \frac{\sin(z\nabla\_2)}{\nabla\_2} \nabla^2 (F + w) \end{cases} \tag{34}$$

Considering the free boundary condition, the shear stress at the plate surface is zero. The equalities can be given from Equations (32) and (33) as follows:

$$\begin{cases} \cos\left(\frac{\hbar}{2}\nabla\_1\right)\frac{\partial w}{\partial x} - \frac{1}{2}\cos\left(\frac{\hbar}{2}\nabla\_2\right) - T\_2^{-2}\sum\_{j=1}^2 (-1)^{j-1}\cos\left(\frac{\hbar}{2}\nabla\_j\right)\nabla^2 \end{cases} \tag{35}$$

$$\begin{cases} \frac{\partial F}{\partial x} + \frac{\partial w}{\partial x} \bigg) - \frac{1}{2}\cos\left(\frac{\hbar}{2}\nabla\_2\right)\frac{\partial f}{\partial y} = 0 \\\\ \cos\left(\frac{\hbar}{2}\nabla\_1\right)\frac{\partial w}{\partial y} - \frac{1}{2}\cos\left(\frac{\hbar}{2}\nabla\_2\right) - T\_2^{-2}\sum\_{j=1}^2 (-1)^{j-1}\cos\left(\frac{\hbar}{2}\nabla\_j\right)\nabla^2 \end{cases} \tag{36}$$

where *h* is the thickness of the plate.

On the basis of the complex variable function theory, Equations (35) and (36) can be regarded as the real part and imaginary part, which consist of a Riemann condition of the analytic function. Now, the non-homogeneous solution of the equation does not influence the solution of the stress state, so there can be:

$$\cos\left(\frac{h}{2}\nabla\_1\right)w - \left[\frac{1}{2}\cos\left(\frac{h}{2}\nabla\_2\right) + T\_2^{-2}\sum\_{j=1}^2 (-1)^{j-1} \cos\left(\frac{h}{2}\nabla\_j\right)\nabla^2\right](F+w) = 0,\tag{37}$$

$$\cos\left(\frac{h}{2}\nabla\_2\right)f = 0.\tag{38}$$

According to the integral function theory, Equation (38) can be expanded in series as:

$$\prod\_{m=1}^{\infty} \left[ 1 - \frac{h^2 \nabla\_2^2}{\left(2m - 1\right)^2 \pi^2} \right] f = 0. \tag{39}$$

By means of truncating the infinite series, the following second-order elastic wave equation is given as:

$$
\nabla^2 f - \left(\frac{\pi^2}{h^2} + \frac{1}{c\_2^2} \frac{\partial^2}{\partial t^2}\right) f = 0. \tag{40}
$$

According to the free boundary condition that the normal stress at the plate surface is zero, the governing equation for the elastic wave of the plate can be derived from Equations (34) and (37) as:

$$T\_2^{-2} \sum\_{j=1}^{2} (-1)^{j-1} \cos(\frac{h}{2} \nabla\_j) \nabla^2 (F + w) + \frac{1}{2} \cos(\frac{h}{2} \nabla\_2)(F + w) - \cos(\frac{h}{2} \nabla\_2)w = 0,\tag{41}$$

$$\begin{cases} T\_2^{-2} \sum\_{j=1 \atop j=1 \atop \sqrt{\frac{1}{2}}}^2 (-1)^{j-1} \frac{\sin(\frac{k}{2} \nabla\_j)}{\nabla\_j} \nabla^2 \nabla^2 (F + w) \ - \frac{1}{2} \left[ \frac{\sin(\frac{k}{2} \nabla\_1)}{\nabla\_1} - \frac{2 \sin(\frac{k}{2} \nabla\_2)}{\nabla\_2} \right] \nabla^2 (F + w) \\\ - \frac{\sin(\frac{k}{2} \nabla\_1)}{\nabla\_1} \nabla^2 w + \frac{1}{2} T\_2^2 \frac{\sin(\frac{k}{2} \nabla\_1)}{\nabla\_1} w = 0 \end{cases} \tag{42}$$

The operator algebraic equation involving the unknown functions *F* and *w* is deduced by Equations (41) and (42) as:

$$
\mathbf{A} \begin{bmatrix} \ F \\ w \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \tag{43}
$$

where the expression of each operator of **Λ** is

$$\begin{split} \Lambda\_{11} &= T\_2^{-2} \sum\_{j=1}^{2} (-1)^{j-1} \cos(\frac{l}{2} \nabla\_j) \nabla^2 + \frac{1}{2} \cos(\frac{l}{2} \nabla\_2), \\ \Lambda\_{12} &= T\_2^{-2} \sum\_{j=1}^{2} (-1)^{j-1} \cos(\frac{l}{2} \nabla\_j) \nabla^2 - \cos(\frac{l}{2} \nabla\_1) + \frac{1}{2} \cos(\frac{l}{2} \nabla\_2), \\ \Lambda\_{21} &= T\_2^{-2} \sum\_{j=1}^{2} (-1)^{j-1} \frac{\sin(\frac{l}{2} \nabla\_j)}{\nabla\_j} \nabla^2 \nabla^2 - \frac{1}{2} \frac{\sin(\frac{l}{2} \nabla\_1)}{\nabla\_1} \nabla^2 + \frac{\sin(\frac{l}{2} \nabla\_2)}{\nabla\_2} \nabla^2, \\ \Lambda\_{22} &= T\_2^{-2} \sum\_{j=1}^{2} (-1)^{j-1} \frac{\sin(\frac{l}{2} \nabla\_j)}{\nabla\_j} \nabla^2 \nabla^2 - \frac{3}{2} \frac{\sin(\frac{l}{2} \nabla\_1)}{\nabla\_1} \nabla^2 + \frac{\sin(\frac{l}{2} \nabla\_2)}{\nabla\_2} \nabla^2 + \frac{1}{2} T\_2^2 \frac{\sin(\frac{l}{2} \nabla\_1)}{\nabla\_1}. \end{split}$$

The determinant of Equation (23) is:

$$\begin{array}{ll} \det(\mathbf{A}) &= T\_2^{-2} \begin{bmatrix} \frac{\sin(\frac{h}{2}\nabla\_1)}{\nabla\_1} \cos(\frac{h}{2}\nabla\_2) - \frac{\sin(\frac{h}{2}\nabla\_2)}{\nabla\_2} \cos(\frac{h}{2}\nabla\_1) \\ - \left[ \frac{\sin(\frac{h}{2}\nabla\_1)}{\nabla\_1} \cos(\frac{h}{2}\nabla\_2) - \frac{\sin(\frac{h}{2}\nabla\_2)}{\nabla\_2} \cos(\frac{h}{2}\nabla\_1) \right] \nabla^2 \\ + \frac{1}{4} T\_2^2 \frac{\sin(\frac{h}{2}\nabla\_1)}{\nabla\_1} \cos(\frac{h}{2}\nabla\_2) \end{bmatrix} \end{array} \tag{44}$$

The fourth-order differential equation involving the lateral displacement can be expressed as:

$$\det(\Lambda)w = 0.\tag{45}$$

After the truncation of the infinite order operator series, the governing equation for the elastic wave of plates can be derived as:

$$
\nabla^2 \nabla^2 w - \frac{3 - 2\kappa}{2(1 - \kappa)} T\_2^2 \nabla^2 w + \frac{3}{1 - \kappa} T\_2^2 \left(\frac{1}{h^2} + \frac{1}{24} T\_1^2 + \frac{1}{8} T\_2^2\right) w = 0,\tag{46}
$$

$$\left[C\nabla^2\nabla^2w - (2-v)DT\_2^2\nabla^2w + \left[CT\_2^2 + \left(\frac{7}{8} - v\right)DT\_2^4\right]w = 0,\tag{47}$$

where *C* = *<sup>E</sup>*M*<sup>h</sup>* <sup>2</sup>(1+*ν*), *<sup>D</sup>* <sup>=</sup> *<sup>E</sup>*M*<sup>h</sup>* <sup>12</sup>(1−*ν*2) , *E*<sup>M</sup> is Young modulus.

Without loss of generality, the solution of the vibration harmonic of the problem is studied. Set:

$$w = \hat{w}\mathbf{e}^{-i\omega t}, \quad F = \hat{F}\mathbf{e}^{-i\omega t}, \; f = \hat{f}\mathbf{e}^{-i\omega t},\tag{48}$$

where *ω* is the angular frequency of plate bending, and i is an imaginary unit.

In the following analysis, the time factor and the symbol '~' in the generalized displacement functions are left out. Taking Equation (48) into Equation (47), the following equations can be expressed as:

$$
\nabla^2 \nabla^2 w + \frac{3 - 2\kappa}{2(1 - \kappa)} k\_2^2 \nabla^2 w + \frac{3}{4(1 - \kappa)} k\_2^2 \left( \frac{\kappa k\_2^2}{6} + \frac{k\_2^2}{2} - \frac{4}{h^2} \right) w = 0,\tag{49}
$$

$$\prod\_{j=1}^{2} \left(\nabla^2 + \mathfrak{a}\_j^2\right) w = 0,\tag{50}$$

where *αj*(*j* = 1, 2) are scattering wave numbers, which satisfy the following expression *<sup>α</sup>*<sup>4</sup> <sup>−</sup> (<sup>2</sup> <sup>−</sup> *<sup>v</sup>*)*k*<sup>2</sup> <sup>2</sup>*α*<sup>2</sup> + *<sup>k</sup>*<sup>2</sup> 2 (7−8*v*)*k*<sup>2</sup> 2 <sup>8</sup> <sup>−</sup> <sup>6</sup>(1−*v*) *h*2 = 0.

The scattering numbers on the basis of Mindlin plate theory are determined by the following expression: *<sup>α</sup>*<sup>4</sup> <sup>−</sup> <sup>12</sup> <sup>π</sup><sup>2</sup> <sup>+</sup> <sup>1</sup>−*<sup>v</sup>* 2 *k*2 <sup>2</sup>*α*<sup>2</sup> <sup>+</sup> <sup>6</sup>(<sup>1</sup> <sup>−</sup> *<sup>v</sup>*)*k*<sup>2</sup> 2 *k*2 2 <sup>π</sup><sup>2</sup> <sup>−</sup> <sup>1</sup> *h*2 = 0, where *kj* = *<sup>ω</sup> cj* , (*j* = 1, 2).

The corresponding generalized displacement potential function is:

$$F = F\_1 + F\_2 = \sum\_{j=1}^{2} (\delta\_j - 1) w\_{j\prime} \tag{51}$$

where *δ<sup>j</sup>* are the ratio coefficients of the displacement potential function, and thus *<sup>δ</sup><sup>j</sup>* <sup>=</sup> <sup>16</sup>+<sup>2</sup> *α*2 *<sup>j</sup>* <sup>−</sup>*κk*<sup>2</sup> 2 *h*2 *h*2 .

$$r\_j = \underset{\underline{\phantom{a}}}{\text{s}} + \overline{\left[ (3-2\kappa)a\_j^2 - k\_j^2 \right]} \,\,^j\underline{\phantom{a}}$$

Comparatively, the ratio coefficients of the displacement potential function by Mindlin plate theory are *<sup>δ</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup> <sup>2</sup> 1−*ν α*2 2*h*2 <sup>π</sup>2−*k*<sup>2</sup> 2*h*<sup>2</sup> , *<sup>δ</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup> <sup>2</sup> 1−*ν α*2 1*h*2 <sup>π</sup>2−*k*<sup>2</sup> 2*h*2 .
