Efficient dynamic analysis of delaminated fractional piezoelectric composites

Abstract

In a variety of applications, such as sensing, vibration control, wind turbine blades, automobile panels, submersible vehicles, and structural health monitoring, smart structures which are made up of laminated composite plates sandwiched between two piezoelectric layers are used. The issue of fractional derivatives for delaminated and multilayered piezoelectric composite plate supported on a viscoelastic foundation’s free vibration response has been tackled in this work. A system of differential quadrature and a perturbation method are used to derive and solve the governing equations. Four form functions, the Lagrange interpolation polynomial, the Cardinal Sine function, the Delta Lagrange, and the Regularized Shannon kernel are used in this study to provide new DQM approaches. A numerical system is suggested and tested under a range of support circumstances to determine its accuracy and efficiency, which relies on prior enabling more accurate predictions critical for the design and deployment of intelligent structural systems.

Introduction

A variety of applications, including sensing, vibration control, wind turbine blades, automobile panels, undersea vehicles, aerospace, and structural health monitoring, use smart hybrid structures, which consist of laminated composite plates sandwiched between two piezoelectric layers. As a result, numerous researchers looked into how a laminated piezoelectric composite plate behaved dynamically. Chen and Lü ¹ investigated laminated plate free vibration analysis using three-dimensional elasticity theory. Khdeir ² looked into laminated plate buckling and free vibration. Matsunaga ³ and Wu and Chen ⁴ used the power series expansion approach to investigate the buckling and vibration of a laminated composite plate. Hachemi ⁵ used the finite element technique (FEM) version p to examine the dynamic behavior of laminated composite sandwich plates with changing stiffness. Using FEM, Kumar and Kumar ⁶ examined the free vibration of laminated composite plates. Wang et al. ⁷ used first order shear deformation theory to study the buckling and free vibration of graphene reinforced composite laminated plates. Additionally, the cylindrical bending vibration of a laminated plate made of uniform layers of piezoelectric composite with varying widths and thicknesses was examined by Vel et al. ⁸ Zhang et al. ⁹ used the differential quadrature method (DQM) to examine the oscillatory behavior of a laminated piezoelectric composite plate. Using DQM and the theories of elasticity and piezoelectricity, Ragb et al. ¹⁰ investigated the oscillatory behavior of piezoelectric composite plates. A piezoelectric composite plate’s oscillatory behavior was precisely measured by Heyliger and Brooks. ¹¹ A piezoelectric laminated composite plate’s static and oscillatory properties were examined by Feri et al. ¹² Piezoelectric laminated materials were examined by Kumari et al. ¹³ More researchers^14–17 used functionally graded materials in their researches. Feri et al. ¹⁸ explored an analytical method for solving three dimensional bending behavior of a functionally graded (FG) material layer covered by piezoelectric layers and subjected to uniform transverse pressure and an electric field. Feri et al. ¹⁹ used state space method and Fourier expansion to investigate thermo viscoelastic response of a sandwich FG plate covered by two piezoelectric layers and it subjected to an electro-thermal load. Ragb and Matbuly ²⁰ studied nonlinear oscillatory behavior of plate resting on Winkler Pasternak foundation using DQM. Zenkour and Alghanmi ²¹ solved a problem of FG plate supported by Winkler–Pasternak foundations attached with an actuator using the Navier method. Sui et al. ²² discussed dynamic behavior of thick plates resting on an elastic foundation, based on the Mindlin plate theory, using FEM. Ye et al. ²³ used FEM to study bending and dynamic analysis of laminated cylindrical panel supported on Pasternak foundation. Chanda et al. ²⁴ studied static and dynamic analysis of laminated piezoelectric composite plate supported on Pasternak foundation. Yas et al. ²⁵ discussed free vibration of FG piezoelectric plates supporting on Pasternak foundation using state space method in the direction of thickness and DQM in the radial direction.

Viscoelastic foundations play a central role in improving the performance and resilience of various engineering structures. Olunloyo et al. ²⁶ proposed an exact solution for vibration problem of a hollow beam conveying fluid on a viscoelastic foundation. Goodarzi et al. ²⁷ discussed vibration of a FG nanoplate supported by visco-Pasternak foundation using modified strain gradient and modified couple stress theories. Rahmani et al. ²⁸ discussed vibration of a laminated composite plate supporting on viscoelastic foundation using DQM.

Delamination is a split between neighboring layers, known as delamination, can be caused by a number of things, including as mechanical stress, heat loading, or environmental deterioration.^29–32 It also spreads along the interface between neighboring layers. One crucial aspect of preserving structural integrity is identifying and fixing delamination. Li and Halim ³³ used the Green function method to investigate the vibration of a delaminated beam. Using FEM, Bhardwaj et al. ³⁴ investigated the vibration and static analysis of a delaminated composite plate. Using FEM, Verenkar et al. ³⁵ investigated the free vibration of a multi-layered delaminated composite plate. Dileep Kumar et al. ³⁶ investigated the dynamic behavior of a delaminated plate using FEM and first order shear deformation theory. FEM was utilized by Keshava Kumar et al. ³⁷ to examine the free vibration of a delaminated composite plate. Moorthy and Marappan ³⁸ used an experimental modality to assess delamination in composite laminates. FEM was utilized by Hirwani et al. ³⁹ to investigate the dynamic analysis of a curved panel made of delaminated composite. FEM was utilized by Ju et al. ⁴⁰ to investigate the dynamic analysis of a delaminated composite. Although many researchers studied nano beams and composite plates on elastic foundations,^41–45 the effect of fractional derivatives have not been investigated.

The Lagrange polynomial interpolation, the Cardinal Sine function, the Delta Lagrange, and the Regularized Shannon kernel are the four shape functions that are used in this paper’s new DQM technique. With few nodal sites, these methods yield precise results.^46–49 Piezoelectricity theory and three-dimensional theory of elasticity is used to create the governing equations supported with fractional derivatives in x-direction. For every differential quadrature approach, MATLAB codes are created. To compare the results with exact and numerical solutions, a comparison analysis is presented. Thus, it is possible to demonstrate the convergence, precision, and effectiveness of the methods that have been provided. Additionally, a parametric research is presented to examine how the frequencies and mode shapes are affected by the length to thickness ratio, damper modulus parameter, delamination length, delamination position, and supporting cases.

Formulation of the problem

Consider a laminated plate consisting of composite and piezoelectric layers which are bounded together. Each layer has a depth of $h_i$ , where 1 ⩽ i ⩽ L, where L is total number of layers. The arrangement of the layers places the 1st layer at the bottom, followed by the 2nd layer above it, and potentially more layers continuing upwards. The overall plate possesses dimensions of length a, width b, and thickness h, $(0\le x\le a,0\le y\le b,0\le z\le h)$ . A delamination has appeared in the x-direction, beginning from a point at a distance of $a_1$ from the left end and continuing to another point at a distance of $a_3$ from the right end edge. In addition, this delamination has a length $a_2$ as shown in Figure 1.

Figure 1.

Laminated plate with a delamination.

Based on the principles of three-dimensional elasticity theory and piezoelectricity theory, the governing equations can be written as^9,50:

$$\frac{\partial {\sigma }_x}{\partial \mathrm{x}}+\frac{\partial {\mathrm{\tau }}_{\mathrm{xy}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{\tau }}_{\mathrm{xz}}}{\partial \mathrm{z}}=\mathrm{\rho }\frac{{\partial }^2\mathrm{u}}{\partial {\mathrm{t}}^2}$$

(1)

$$\frac{\partial {\mathrm{\tau }}_{\mathrm{xy}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{\sigma }}_{\mathrm{y}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{\tau }}_{\mathrm{yz}}}{\partial \mathrm{z}}=\mathrm{\rho }\frac{{\partial }^2\mathrm{v}}{\partial {\mathrm{t}}^2}$$

(2)

$$\frac{\partial {\mathrm{\tau }}_{\mathrm{x}\mathrm{z}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{\tau }}_{\mathrm{y}\mathrm{z}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{\sigma }}_{\mathrm{z}}}{\partial \mathrm{z}}+{\mathrm{C}}_{\mathrm{d}}\frac{\partial \mathrm{w}}{\partial \mathrm{t}}=\mathrm{\rho }\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{t}}^2}$$

(3)

$$\frac{\partial {\mathrm{D}}_{\mathrm{x}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{D}}_{\mathrm{y}}}{\partial \mathrm{y}}+\frac{\partial {\mathrm{D}}_{\mathrm{z}}}{\partial \mathrm{z}}=0$$

(4)

In these equations, (x, y, z) represent the reference coordinate system. In addition (u, v, w) refer to the displacements in x, y, z directions respectively, t represents time, $\mathrm{\rho }$ denotes mass density, $({\mathrm{\sigma }}_{\mathrm{x}},{\mathrm{\sigma }}_{\mathrm{y}},{\mathrm{\sigma }}_{\mathrm{z}})$ represent stresses, $({\mathrm{\tau }}_{\mathrm{xy}},{\mathrm{\tau }}_{\mathrm{xz}},{\mathrm{\tau }}_{\mathrm{yz}})$ refer to shear stresses, ${\mathrm{C}}_{\mathrm{d}}$ is damper modulus parameter, and $({\mathrm{D}}_{\mathrm{x}},{\mathrm{D}}_{\mathrm{y}},{\mathrm{D}}_{\mathrm{z}})$ are the electric displacements.

The formulations which define the constitutive properties for a cross-ply laminate of an orthotropic material can be expressed as^9,50:

$$\left[\begin{array}{c}\begin{array}{c}{\mathrm{\sigma }}_{\mathrm{x}}\\ {\mathrm{\sigma }}_{\mathrm{y}}\\ {\mathrm{\sigma }}_{\mathrm{z}}\end{array}\\ \begin{array}{c}{\mathrm{\tau }}_{\mathrm{yz}}\\ {\mathrm{\tau }}_{\mathrm{xz}}\\ {\mathrm{\tau }}_{\mathrm{xy}}\end{array}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{ccc}{\mathrm{c}}_{11}& {\mathrm{c}}_{12}& {\mathrm{c}}_{13}\\ {\mathrm{c}}_{12}& {\mathrm{c}}_{22}& {\mathrm{c}}_{23}\\ {\mathrm{c}}_{13}& {\mathrm{c}}_{23}& {\mathrm{c}}_{33}\end{array}& \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}{\mathrm{c}}_{44}& 0& 0\\ 0& {\mathrm{c}}_{55}& 0\\ 0& 0& {\mathrm{c}}_{66}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\mathrm{u}}_{\mathrm{x}}\\ {\mathrm{v}}_{\mathrm{y}}\\ {\mathrm{w}}_{\mathrm{z}}\end{array}\\ \begin{array}{c}{\mathrm{w}}_{\mathrm{y}}+{\mathrm{v}}_{\mathrm{z}}\\ {\mathrm{w}}_{\mathrm{x}}+{\mathrm{u}}_{\mathrm{z}}\\ {\mathrm{v}}_{\mathrm{x}}+{\mathrm{u}}_{\mathrm{y}}\end{array}\end{array}\right]+\left[\begin{array}{c}\begin{array}{ccc}0& 0& {\mathrm{e}}_1\\ 0& 0& {\mathrm{e}}_2\\ 0& 0& {\mathrm{e}}_3\end{array}\\ \begin{array}{ccc}0& {\mathrm{e}}_4& 0\\ {\mathrm{e}}_5& 0& 0\\ 0& 0& 0\end{array}\end{array}\right]\left[\begin{array}{c}{\mathrm{\phi }}_{\mathrm{x}}\\ {\mathrm{\phi }}_{\mathrm{y}}\\ {\mathrm{\phi }}_{\mathrm{z}}\end{array}\right]$$

(5)

$$\left[\begin{array}{c}{\mathrm{D}}_{\mathrm{x}}\\ {\mathrm{D}}_{\mathrm{y}}\\ {\mathrm{D}}_{\mathrm{z}}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ {\mathrm{e}}_1& {\mathrm{e}}_2& {\mathrm{e}}_3\end{array}& \begin{array}{ccc}0& {\mathrm{e}}_5& 0\\ {\mathrm{e}}_4& 0& 0\\ 0& 0& 0\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\mathrm{u}}_{\mathrm{x}}\\ {\mathrm{v}}_{\mathrm{y}}\\ {\mathrm{w}}_{\mathrm{z}}\end{array}\\ \begin{array}{c}{\mathrm{w}}_{\mathrm{y}}+{\mathrm{v}}_{\mathrm{z}}\\ {\mathrm{w}}_{\mathrm{x}}+{\mathrm{u}}_{\mathrm{z}}\\ {\mathrm{v}}_{\mathrm{x}}+{\mathrm{u}}_{\mathrm{y}}\end{array}\end{array}\right]-\left[\begin{array}{ccc}{\mathrm{\eta }}_1& 0& 0\\ 0& {\mathrm{\eta }}_2& 0\\ 0& 0& {\mathrm{\eta }}_3\end{array}\right]\left[\begin{array}{c}{\mathrm{\phi }}_{\mathrm{x}}\\ {\mathrm{\phi }}_{\mathrm{y}}\\ {\mathrm{\phi }}_{\mathrm{z}}\end{array}\right]$$

(6)

where ${\mathrm{c}}_{11},{\mathrm{c}}_{12},\dots \dots$ can be written as follows^1,3:

$$\left\{\begin{array}{c}{c}_{11}\\ c_{22}\\ c_{12}\\ c_{13}\\ c_{23}\\ c_{33}\\ c_{44}\\ c_{55}\\ c_{66}\end{array}\right\}=\left\{\begin{array}{c}\overline{c_{11}}{\cos }^4\mathrm{\theta }+\overline{c_{22}}{\sin }^4\mathrm{\theta }\\ \overline{c_{11}}{\sin }^4\mathrm{\theta }+\overline{c_{22}}{\cos }^4\mathrm{\theta }\\ \overline{c_{12}}\left({\cos }^4\mathrm{\theta }+{\sin }^4\mathrm{\theta }\right)\\ \overline{c_{13}}{\cos }^2\mathrm{\theta }+\overline{c_{23}}{\sin }^2\mathrm{\theta }\\ \overline{c_{13}}{\sin }^2\mathrm{\theta }+\overline{c_{23}}{\cos }^2\mathrm{\theta }\\ \overline{c_{33}}\\ \overline{c_{44}}{\cos }^2\mathrm{\theta }+\overline{c_{55}}{\sin }^2\mathrm{\theta }\\ \overline{c_{44}}{\sin }^2\mathrm{\theta }+\overline{c_{55}}{\cos }^2\mathrm{\theta }\\ \overline{c_{66}}\end{array}\right\},\left\{\begin{array}{c}\overline{c_{11}}\\ \overline{c_{12}}=\overline{c_{21}}\\ \overline{c_{13}}=\overline{c_{31}}\\ \overline{c_{22}}\\ \overline{c_{23}}=\overline{c_{32}}\\ \overline{c_{33}}\\ \overline{c_{44}}\\ \overline{c_{55}}\\ \overline{c_{66}}\end{array}\right\}=\left\{\begin{array}{c}\frac{E_1\left(E_2-{\mathrm{\nu }}_{23}^2{E}_3\right)}{\mu E_2}\\ \frac{{\mathrm{\nu }}_{12}{E}_2+{\mathrm{\nu }}_{13}{\mathrm{\nu }}_{23}{E}_3}{\mu }\\ \frac{E_3\left({\mathrm{\nu }}_{13}-{\mathrm{\nu }}_{12}{\mathrm{\nu }}_{23}\right)}{\mu }\\ \frac{E_2\left(E_1-{\mathrm{\nu }}_{13}^2{E}_3\right)}{\mu E_1}\\ \frac{E_3\left({\mathrm{\nu }}_{23}{E}_1-{\mathrm{\nu }}_{12}{\mathrm{\nu }}_{13}{E}_2\right)}{\mu E_1}\\ \frac{E_3\left(E_1-{\mathrm{\nu }}_{12}^2{E}_2\right)}{\mu E_1}\\ {\mathrm{G}}_{23}\\ {\mathrm{G}}_{13}\\ {\mathrm{G}}_{12}\end{array}\right\}$$

(7)

$$\mu =1-{\mathrm{\nu }}_{12}^2\frac{E_2}{E_1}-{\mathrm{\nu }}_{23}^2\frac{E_3}{E_2}-{\mathrm{\nu }}_{13}^2\frac{E_3}{E_1}-2{\mathrm{\nu }}_{12}{\mathrm{\nu }}_{23}{\mathrm{\nu }}_{13}\frac{E_3}{E_1}$$

(8)

In these relations, $E_{\mathrm{i}},G_{\mathrm{i}\mathrm{j}},{\mathrm{\nu }}_{\mathrm{i}\mathrm{j}}$ denote to Young’s moduli, shear moduli and Poisson’s ratios. $\mathrm{\phi },{\mathrm{e}}_{\mathrm{i}},{\mathrm{\eta }}_{\mathrm{i}}$ refer to electrical potential, effective piezoelectric constant, and effective dielectric constant respectively. The angle θ represents the angle between the x-direction and the local axis 1 of a layer.

Assuming harmonic behavior in the free vibration of the plate, the field quantities can be written as:

$$\mathrm{f}(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t})=\mathrm{F}{\mathrm{e}}^{\mathrm{i}\mathrm{\omega }\mathrm{t}},(\mathrm{i}=\sqrt{-1})$$

(9)

Here $\mathrm{\omega }$ represents natural frequency of the plate. F is amplitude for f.

Hence, we apply the equation (9) into the equations (1) to (6), one can write the governing equations as:

$$\begin{array}{l}{\mathrm{c}}_{11}\frac{{\partial }^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{U}}{\partial {\mathrm{x}}^{\mathrm{\alpha },\mathrm{\epsilon }}}+{\mathrm{c}}_{66}\frac{{\partial }^2\mathrm{U}}{\partial {\mathrm{y}}^2}+{\mathrm{c}}_{55}\frac{{\partial }^{\mathrm{\beta }}\mathrm{U}}{\partial {\mathrm{z}}^{\mathrm{\beta }}}+({\mathrm{c}}_{12}+{\mathrm{c}}_{66})\frac{{\partial }^2\mathrm{V}}{\partial \mathrm{x}\mathrm{y}}\\ +({\mathrm{c}}_{13}+{\mathrm{c}}_{55})\frac{{\partial }^2\mathrm{W}}{\partial \mathrm{x}\mathrm{z}}+({\mathrm{e}}_1+{\mathrm{e}}_5)\frac{{\partial }^2\mathrm{\varphi }}{\partial \mathrm{x}\mathrm{z}}=-\mathrm{\rho }{\mathrm{\omega }}^2\mathrm{U}\end{array}$$

(10)

$$\begin{array}{l}({\mathrm{c}}_{12}+{\mathrm{c}}_{66})\frac{{\partial }^2\mathrm{U}}{\partial \mathrm{xy}}+{\mathrm{c}}_{66}\frac{{\partial }^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{V}}{\partial {\mathrm{x}}^{\mathrm{\alpha },\mathrm{\epsilon }}}+{\mathrm{c}}_{22}\frac{{\partial }^2\mathrm{V}}{\partial {\mathrm{y}}^2}\\ +{\mathrm{c}}_{44}\frac{{\partial }^{\mathrm{\beta }}\mathrm{V}}{\partial {\mathrm{z}}^{\mathrm{\beta }}}+({\mathrm{c}}_{23}+{\mathrm{c}}_{44})\frac{{\partial }^2\mathrm{W}}{\partial \mathrm{yz}}\\ +({\mathrm{e}}_2+{\mathrm{e}}_4)\frac{{\partial }^2\mathrm{\varphi }}{\partial \mathrm{yz}}=-\mathrm{\rho }{\mathrm{\omega }}^2\mathrm{V}\end{array}$$

(11)

$$\begin{array}{l}({\mathrm{c}}_{13}+{\mathrm{c}}_{55})\frac{{\partial }^2\mathrm{U}}{\partial \mathrm{x}\mathrm{z}}+({\mathrm{c}}_{23}+{\mathrm{c}}_{44})\frac{{\partial }^2\mathrm{V}}{\partial \mathrm{y}\mathrm{z}}+({\mathrm{c}}_{55})\frac{{\partial }^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{W}}{\partial {\mathrm{x}}^{\mathrm{\alpha },\mathrm{\epsilon }}}\\ +({\mathrm{c}}_{44})\frac{{\partial }^2\mathrm{W}}{\partial {\mathrm{y}}^2}+{\mathrm{c}}_{33}\frac{{\partial }^{\mathrm{\beta }}\mathrm{W}}{\partial {\mathrm{z}}^{\mathrm{\beta }}}+{\mathrm{e}}_5\frac{{\partial }^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{\varphi }}{\partial {\mathrm{x}}^{\mathrm{\alpha },\mathrm{\epsilon }}}+{\mathrm{e}}_4\frac{{\partial }^2\mathrm{\varphi }}{\partial {\mathrm{y}}^2}\\ +{\mathrm{e}}_3\frac{{\partial }^{\mathrm{\beta }}\mathrm{\varphi }}{\partial {\mathrm{z}}^{\mathrm{\beta }}}+{\mathrm{C}}_{\mathrm{d}}(\mathrm{i}\mathrm{\omega })\mathrm{W}=-\mathrm{\rho }{\mathrm{\omega }}^2\mathrm{W}\end{array}$$

(12)

$$\begin{array}{l}({\mathrm{e}}_1+{\mathrm{e}}_5)\frac{{\partial }^2\mathrm{U}}{\partial \mathrm{x}\mathrm{z}}+({\mathrm{e}}_2+{\mathrm{e}}_4)\frac{{\partial }^2\mathrm{V}}{\partial \mathrm{y}\mathrm{z}}+{\mathrm{e}}_5\frac{{\partial }^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{W}}{\partial {\mathrm{x}}^{\mathrm{\alpha },\mathrm{\epsilon }}}+{\mathrm{e}}_4\frac{{\partial }^2\mathrm{W}}{\partial {\mathrm{y}}^2}\\ +{\mathrm{e}}_3\frac{{\partial }^{\mathrm{\beta }}\mathrm{W}}{\partial {\mathrm{z}}^{\mathrm{\beta }}}-{\mathrm{\eta }}_1\frac{{\partial }^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{\varphi }}{\partial {\mathrm{x}}^{\mathrm{\alpha },\mathrm{\epsilon }}}-{\mathrm{\eta }}_2\frac{{\partial }^2\mathrm{\varphi }}{\partial {\mathrm{y}}^2}-{\mathrm{\eta }}_3\frac{{\partial }^{\mathrm{\beta }}\mathrm{\varphi }}{\partial {\mathrm{z}}^{\mathrm{\beta }}}=0\end{array}$$

(13)

The boundary conditions can be written as ⁹ :

For a simply supported end:

$$\mathrm{W}=\mathrm{V}={\mathrm{\sigma }}_{\mathrm{x}}=\mathrm{\varphi }=0,\mathrm{at}\mathrm{x}=0,\mathrm{a}$$

(14)

$$\mathrm{W}=\mathrm{U}={\mathrm{\sigma }}_{\mathrm{y}}=\mathrm{\varphi }=0,\mathrm{at}\mathrm{y}=0,\mathrm{b}$$

(15)

For a clamped supported end:

$$\mathrm{U}=\mathrm{V}=\mathrm{W}=\mathrm{\varphi }=0,\mathrm{at}\mathrm{x}=0,\mathrm{a},\mathrm{y}=0,\mathrm{b}$$

(16)

For a free supported end:

$${\mathrm{\sigma }}_{\mathrm{x}}={\mathrm{\tau }}_{\mathrm{x}\mathrm{z}}={\mathrm{\tau }}_{\mathrm{x}\mathrm{y}}=\mathrm{\varphi }=0,\mathrm{at}\mathrm{x}=0,\mathrm{a}$$

(17)

$${\mathrm{\sigma }}_{\mathrm{y}}={\mathrm{\tau }}_{\mathrm{y}\mathrm{z}}={\mathrm{\tau }}_{\mathrm{x}\mathrm{y}}=\mathrm{\varphi }=0,\mathrm{at}\mathrm{y}=0,\mathrm{b}$$

(18)

For upper and lower surfaces of the plate, boundary condition can be written as:

For closed circuit (grounded):

$${\mathrm{\sigma }}_{\mathrm{z}}={\mathrm{\tau }}_{\mathrm{x}\mathrm{z}}={\mathrm{\tau }}_{\mathrm{y}\mathrm{z}}=\mathrm{\varphi }=0,\mathrm{a}\mathrm{t}\mathrm{z}=0,\mathrm{h}$$

(19)

For open circuit:

$${\mathrm{\sigma }}_{\mathrm{z}}={\mathrm{\tau }}_{\mathrm{x}\mathrm{z}}={\mathrm{\tau }}_{\mathrm{y}\mathrm{z}}={\mathrm{D}}_{\mathrm{z}}=0,\mathrm{a}\mathrm{t}\mathrm{z}=0,\mathrm{h}$$

(20)

At the interface locations, the continuity conditions can be expressed as:

$$\begin{array}{l}{(\mathrm{U})}_{\mathrm{p}}={(\mathrm{U})}_{\mathrm{p}+1},{(\mathrm{V})}_{\mathrm{p}}={(\mathrm{V})}_{\mathrm{p}+1},\\ {(\mathrm{W})}_{\mathrm{p}}={(\mathrm{W})}_{\mathrm{p}+1},{(\mathrm{\varphi })}_{\mathrm{p}}={(\mathrm{\varphi })}_{\mathrm{p}+1}\\ {({\mathrm{\sigma }}_{\mathrm{z}})}_{\mathrm{p}}={({\mathrm{\sigma }}_{\mathrm{z}})}_{\mathrm{p}+1},{({\mathrm{\tau }}_{\mathrm{x}\mathrm{z}})}_{\mathrm{p}}={({\mathrm{\tau }}_{\mathrm{x}\mathrm{z}})}_{\mathrm{p}+1},\\ {({\mathrm{\tau }}_{\mathrm{y}\mathrm{z}})}_{\mathrm{p}}={({\mathrm{\tau }}_{\mathrm{y}\mathrm{z}})}_{\mathrm{p}+1},{({\mathrm{D}}_{\mathrm{z}})}_{\mathrm{p}}={({\mathrm{D}}_{\mathrm{z}})}_{\mathrm{p}+1}\end{array}$$

(21)

Here, p = (1, L−1), L denotes to the number of layers.

Method of solution

Polynomial DQM (classical DQM)

In this approach, the Lagrange interpolation polynomial is used as the shape function. The nodal points can be uniformly or non-uniformly distributed. To get the non-uniform distribution of the points, one may use the Chebyshev Gauss Lobatto discretization as follows ⁵¹ :

$${\mathrm{x}}_{\mathrm{i}}={\mathrm{x}}_1+\frac{{\mathrm{x}}_{\mathrm{N}}-{\mathrm{x}}_1}{2}\left(1-\cos \left(\frac{\mathrm{i}-1}{\mathrm{N}-1}\mathrm{\pi }\right)\right),\mathrm{i}=1,2,··,\mathrm{N}$$

(22)

This technique can approximate the field quantities as follows^51–55:

$$\mathrm{f}({\mathrm{x}}_{\mathrm{i}})=\sum _{\mathrm{j}=1}^{\mathrm{N}}\frac{1}{{\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}}\left(\frac{{\prod }_{\mathrm{k}=1,\mathrm{k}\ne \mathrm{i}}^{\mathrm{N}}({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{k}})}{{\prod }_{\mathrm{k}=1,\mathrm{k}\ne \mathrm{j}}^{\mathrm{N}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}})}\right)\mathrm{f}({\mathrm{x}}_{\mathrm{j}}),(\mathrm{i}=1,\mathrm{N})$$

(23)

For the 1st derivative, the weighting coefficients can be expressed as^51–56:

$${\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1=\{\begin{array}{c}\frac{1}{{\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}}\left(\frac{{\mathrm{P}}_{\mathrm{i}}}{{\mathrm{P}}_{\mathrm{j}}}\right),\mathrm{i}\ne \mathrm{j}\\ -\sum _{\mathrm{j}=1,\mathrm{j}\ne \mathrm{i}}^{\mathrm{N}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1,\mathrm{i}=\mathrm{j}\end{array}$$

(24)

Here, ${\mathrm{P}}_{\mathrm{i}}$ is defined as:

$${\mathrm{P}}_{\mathrm{i}}=\prod _{\mathrm{k}=1,\mathrm{k}\ne \mathrm{i}}^{\mathrm{N}}({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{k}}),\mathrm{k}=1,2,··,\mathrm{N}$$

(25)

For the nth derivative, the weighting coefficients may be represented as ⁵⁵ :

$$[{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{n}}]=[{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1][{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{n}-1}]$$

(26)

Sinc DQM

The Sinc DQM uses the Cardinal Sine function as a shape function. The points are uniformly distributed as uniform points with number N. The field quantities and their derivatives are transformed to weighted linear sums of the nodal values f_i,(i = −N : N), as follows^10,55:

$$\begin{array}{l}\mathrm{f}({\mathrm{x}}_{\mathrm{i}})=\sum _{\mathrm{j}=-\mathrm{N}}^{\mathrm{N}}\frac{\sin [\mathrm{\pi }({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}})/{\mathrm{h}}_{\mathrm{m}}]}{\mathrm{\pi }({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}})/{\mathrm{h}}_{\mathrm{m}}}\mathrm{f}({\mathrm{x}}_{\mathrm{j}}),\\ (\mathrm{i}=-\mathrm{N},\mathrm{N}),{\mathrm{h}}_{\mathrm{m}}>0\end{array}$$

(27)

Weighting coefficients for the 1st and the 2nd derivatives are represented as follows^10,55:

$$\begin{array}{l}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1=\{\begin{array}{c}\frac{{(-1)}^{\mathrm{i}-\mathrm{j}}}{{\mathrm{h}}_{\mathrm{m}}(\mathrm{i}-\mathrm{j})}\mathrm{i}\ne \mathrm{j}\\ 0\mathrm{i}=\mathrm{j}\end{array},\\ {\mathrm{C}}_{\mathrm{i}\mathrm{j}}^2=\{\begin{array}{c}\frac{2{(-1)}^{\mathrm{i}-\mathrm{j}+1}}{{\mathrm{h}}_{\mathrm{m}}^2{(\mathrm{i}-\mathrm{j})}^2}\mathrm{i}\ne \mathrm{j}\\ -\frac{{\mathrm{\pi }}^2}{3{\mathrm{h}}_{\mathrm{m}}^2}\mathrm{i}=\mathrm{j}\end{array}\end{array}$$

(28)

Here, ${\mathrm{h}}_{\mathrm{m}}$ refers to the interval between two points, $({\mathrm{h}}_{\mathrm{x}}={\mathrm{x}}_{\mathrm{i}+1}-{\mathrm{x}}_{\mathrm{i}}),({\mathrm{h}}_{\mathrm{y}}={\mathrm{y}}_{\mathrm{i}+1}-{\mathrm{y}}_{\mathrm{i}}),({\mathrm{h}}_{\mathrm{z}}={\mathrm{z}}_{\mathrm{i}+1}-{\mathrm{z}}_{\mathrm{i}}).$

Discrete singular convolution DQM

The singular convolution can be written as follows^10,55:

$$\mathrm{F}(\mathrm{t})=(\mathrm{T}*\mathrm{\zeta })(\mathrm{t})={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{T}(\mathrm{t}-\mathrm{x})\mathrm{\zeta }(\mathrm{x})\mathrm{d}\mathrm{x}$$

(29)

Where, $\mathrm{T}(\mathrm{t}-\mathrm{x})$ refers to a singular kernel.

This technique depends on many kernels used as shape functions. The points are uniformly distributed as uniform points with number N. This technique uses narrow bandwidth (2M + 1) to can control the accuracy of the results. There are two common kernels known as Delta Lagrange Kernel (DLK) and Regularized Shannon kernel (RSK).

Delta Lagrange Kernel (DLK)

Delta Lagrange Kernel (DLK) may be used as a shape function. So, the field quantities can be written as:

$$\begin{array}{l}\mathrm{f}({\mathrm{x}}_{\mathrm{i}})=\sum _{\mathrm{j}=-\mathrm{M}}^{\mathrm{M}}\frac{1}{({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}})}\left(\frac{{\prod }_{\mathrm{k}=-\mathrm{M},\mathrm{k}\ne \mathrm{i}}^{\mathrm{M}}({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{k}})}{{\prod }_{\mathrm{k}=-\mathrm{M},\mathrm{k}\ne \mathrm{j}}^{\mathrm{M}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{k}})}\right)\mathrm{f}({\mathrm{x}}_{\mathrm{j}}),\\ (\mathrm{i}=-\mathrm{N},\mathrm{N}),\mathrm{M}\ge 1\end{array}$$

(30)

The weighting coefficients of the 1st and the 2nd derivatives can be defined as^10,55:

$$\begin{array}{l}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1=\{\begin{array}{c}\frac{1}{{\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}}\left(\frac{{\mathrm{P}}_{\mathrm{i}}}{{\mathrm{P}}_{\mathrm{j}}}\right),\mathrm{i}\ne \mathrm{j}\\ -\sum _{\mathrm{j}=-\mathrm{M},\mathrm{j}\ne \mathrm{i}}^{\mathrm{M}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1,\mathrm{i}=\mathrm{j}\end{array},\\ {\mathrm{C}}_{\mathrm{i}\mathrm{j}}^2=\{\begin{array}{c}2{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1\left({\mathrm{C}}_{\mathrm{i}\mathrm{i}}^1-\frac{1}{{\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}}\right)\mathrm{i}\ne \mathrm{j}\\ -\sum _{\mathrm{j}=-\mathrm{M},\mathrm{j}\ne \mathrm{i}}^{\mathrm{M}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^2\mathrm{i}=\mathrm{j}\end{array}\end{array}$$

(31)

where ${\mathrm{P}}_{\mathrm{i}}$ is defined as:

$${\mathrm{P}}_{\mathrm{i}}=\prod _{\mathrm{k}=-\mathrm{M},\mathrm{k}\ne \mathrm{i}}^{\mathrm{M}}({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{k}})$$

(32)

Regularized Shannon kernel (RSK)

Regularized Shannon kernel (RSK) is used as a shape function as:

$$\begin{array}{l}\mathrm{f}({\mathrm{x}}_{\mathrm{i}})=\sum _{\mathrm{j}=-\mathrm{M}}^{\mathrm{M}}\left(\frac{\sin [\mathrm{\pi }({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}})/{\mathrm{h}}_{\mathrm{m}}]}{\mathrm{\pi }({\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}})/{\mathrm{h}}_{\mathrm{m}}}{\mathrm{e}}^{-\left(\frac{{{\mathrm{h}}_{\mathrm{m}}}^2{(\mathrm{i}-\mathrm{j})}^2}{2{\mathrm{\sigma }}^2}\right)}\right)\mathrm{f}({\mathrm{x}}_{\mathrm{j}}),\\ (\mathrm{i}=-\mathrm{N},\mathrm{N}),\mathrm{\sigma }=(\mathrm{r}*{\mathrm{h}}_{\mathrm{m}})>0\end{array}$$

(33)

The 1st and the 2nd derivatives have weighting coefficients as follows^10,53,56

$$\begin{array}{l}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1=\{\begin{array}{c}\frac{{(-1)}^{\mathrm{i}-\mathrm{j}}}{{\mathrm{h}}_{\mathrm{m}}(\mathrm{i}-\mathrm{j})}{\mathrm{e}}^{-\left(\frac{{\mathrm{h}}_{\mathrm{m}}^2{(\mathrm{i}-\mathrm{j})}^2}{2{\mathrm{\sigma }}^2}\right)}\mathrm{i}\ne \mathrm{j}\\ 0\mathrm{i}=\mathrm{j}\end{array},\\ {\mathrm{C}}_{\mathrm{i}\mathrm{j}}^2=\{\begin{array}{c}\left(\frac{2{(-1)}^{\mathrm{i}-\mathrm{j}+1}}{{\mathrm{h}}_{\mathrm{m}}^2{(\mathrm{i}-\mathrm{j})}^2}+\frac{1}{{\mathrm{\sigma }}^2}\right){\mathrm{e}}^{-\left(\frac{{\mathrm{h}}_{\mathrm{m}}^2{(\mathrm{i}-\mathrm{j})}^2}{2{\mathrm{\sigma }}^2}\right)}\mathrm{i}\ne \mathrm{j}\\ -\frac{1}{{\mathrm{\sigma }}^2}-\frac{{\mathrm{\pi }}^2}{3{\mathrm{h}}_{\mathrm{m}}^2}\mathrm{i}=\mathrm{j}\end{array}\end{array}$$

(34)

where, $\mathrm{\sigma }=\mathrm{r}*{\mathrm{h}}_{\mathrm{m}}>0$ , is the regularization parameter and r is the computational parameter.

Fractional derivatives

Fractional derivative is an important subject in many engineering applications. There are different definitions of fractional derivative, such as Riemann–Liouville, Caputo, and Weyl.

Caputo’s fractional derivative

By assuming (n−1) < α < n, the fractional derivative of f(x) of order α based can be written as ⁴⁶ :

$${\mathrm{D}}_{\mathrm{c}}^{\mathrm{\alpha }}\mathrm{f}(\mathrm{X})=\{\begin{array}{c}\begin{array}{l}\frac{1}{\mathrm{\Gamma }(\mathrm{n}-\mathrm{\alpha })}\underset{\mathrm{c}}{\overset{\mathrm{X}}{\int }}{(\mathrm{X}-\mathrm{x})}^{\mathrm{n}-\mathrm{\alpha }-1}{\mathrm{f}}^{\mathrm{k}}(\mathrm{x})\mathrm{d}\mathrm{x},\\ (\mathrm{n}-1)<\mathrm{\alpha }<\mathrm{n}\end{array}\\ \frac{{\mathrm{d}}^{\mathrm{n}}\mathrm{f}}{\mathrm{d}{\mathrm{x}}^{\mathrm{n}}}\mathrm{\alpha }=\mathrm{n}\end{array}$$

(35)

where c is the lower limit of integration.

Caputo’s fractional derivative can be applied with classical PDQM, Sinc DQM, DLK DQM and RSK DQM by using equation (35) and equation (22) to determine the weighting coefficients ${\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{\alpha }}$ for α ∈ (0, 1] as ⁴⁶ :

$${\mathrm{D}}^{\mathrm{\alpha }}\mathrm{f}(\mathrm{X})=\{\begin{array}{c}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{x}}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{\alpha }}\mathrm{f}({\mathrm{X}}_{\mathrm{j}},\mathrm{x}),0<\mathrm{\alpha }<1\\ {\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{x}}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{x}}\mathrm{f}({\mathrm{X}}_{\mathrm{j}},\mathrm{x})\mathrm{\alpha }=1\end{array}$$

(36)

where

${\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{\alpha }}={\mathrm{A}}^{1-\mathrm{\alpha }}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{x}}-\frac{{\mathrm{C}}_{1,\mathrm{j}}^{\mathrm{x}}}{\mathrm{\Gamma }(2-\mathrm{\alpha })}{(\mathrm{X}-\mathrm{c})}^{1-\mathrm{\alpha }}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{A}}_{\mathrm{i}\mathrm{j}}={\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{x}}-{\mathrm{C}}_{1,\mathrm{j}}^{\mathrm{x}}$

Generalized Caputo’s fractional derivative

Fractional differential operators due to their non-local characteristics, which can take various forms, most effectively describe the system influenced by memory. Therefore, the researchers extend the fractional operators to hold the hidden aspects of the real nonlocal phenomena. The definition of the generalized Caputo fractional derivative of a certain order and operator is presented by^44,45:

$$\begin{array}{l}{\mathrm{D}}_{{\mathrm{c}}_{+}}^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{f}(\mathrm{X})=\frac{{\mathrm{\rho }}^{\mathrm{\alpha }-\mathrm{n}+1}}{\mathrm{\Gamma }(\mathrm{n}-\mathrm{\alpha })}\underset{\mathrm{c}}{\overset{\mathrm{X}}{\int }}{\mathrm{x}}^{\mathrm{\epsilon }-1}{({\mathrm{X}}^{\mathrm{\epsilon }}-{\mathrm{x}}^{\mathrm{\epsilon }})}^{\mathrm{n}-\mathrm{\alpha }-1}{\left({\mathrm{x}}^{1-\mathrm{\epsilon }}\frac{\mathrm{d}}{\mathrm{d}\mathrm{x}}\right)}^{\mathrm{n}}\\ \mathrm{f}(\mathrm{x})\mathrm{d}\mathrm{x},(\mathrm{n}-1)<\mathrm{\alpha }<\mathrm{n},\mathrm{\epsilon }>0,\mathrm{c}\ge 0\end{array}$$

(37)

Henceforth, the solution of equation (37) can be written as:

$${\mathrm{D}}_{{\mathrm{c}}_{+}}^{\mathrm{\alpha },\mathrm{\epsilon }}{({\mathrm{X}}^{\mathrm{\epsilon }}-{\mathrm{x}}^{\mathrm{\epsilon }})}^{\mathrm{n}}={\mathrm{\epsilon }}^{\mathrm{\alpha }}\frac{\mathrm{\Gamma }(\mathrm{n}+1)}{\mathrm{\Gamma }(\mathrm{n}-\mathrm{\alpha }+1)}{({\mathrm{X}}^{\mathrm{\epsilon }}-{\mathrm{c}}^{\mathrm{\epsilon }})}^{\mathrm{n}-\mathrm{\alpha }}$$

(38)

By substituting the generalized Caputo’s fractional derivative in equation (37) into equation (22) of differential quadrature method based on four shape functions to calculate the weighting coefficients ${\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{\alpha },\mathrm{\epsilon }}$ for α ∈ (0, 1] and ε > 0, as ⁵⁵ :

$${\mathrm{D}}_{{\mathrm{c}}_{+}}^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{f}(\mathrm{X})=\{\begin{array}{c}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{x}}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{\alpha },\mathrm{\epsilon }}\mathrm{f}({\mathrm{X}}_{\mathrm{j}},\mathrm{x}),0<\mathrm{\alpha }\le 1,\mathrm{\epsilon }>0\\ {\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{x}}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{(1)}\mathrm{f}({\mathrm{X}}_{\mathrm{j}},\mathrm{x})\mathrm{\alpha }=\mathrm{\epsilon }=1\end{array}$$

(39)

So, the weighting coefficient ${\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{\alpha },\mathrm{\epsilon }}$ is written as ⁵⁶ :

$${\mathrm{C}}_{\mathrm{i}\mathrm{j}}^{\mathrm{\alpha },\mathrm{\epsilon }}={\mathrm{A}}^{1-\mathrm{\alpha }}{\mathrm{\epsilon }}^{\mathrm{\alpha }}{\mathrm{C}}_{\mathrm{i}\mathrm{j}}^1+\frac{{\mathrm{\epsilon }}^{\mathrm{\alpha }}{\mathrm{C}}_{1,\mathrm{j}}^1}{\mathrm{\Gamma }(2-\mathrm{\alpha })}{({\mathrm{X}}^{\mathrm{\epsilon }}-{\mathrm{c}}^{\mathrm{\epsilon }})}^{1-\mathrm{\alpha }}$$

(40)

Discretization of the problem

In this section, we substitute DQM equations with four shape functions as appearing in equations (22) to (40) on equations (10) to (21). Consequently, the problem is transformed to the following non-linear algebraic system:

$$\begin{array}{l}{\mathrm{c}}_{11}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^{\mathrm{\alpha },\mathrm{\epsilon }}}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{66}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^2}^{\mathrm{y}}{\mathrm{U}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{55}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^{\mathrm{\beta }}}^{\mathrm{z}}{\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+({\mathrm{c}}_{12}+{\mathrm{c}}_{66})\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{l}\mathrm{m}\mathrm{k}}\\ +({\mathrm{c}}_{13}+{\mathrm{c}}_{55})\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{l}\mathrm{j}\mathrm{n}}+({\mathrm{e}}_1+{\mathrm{e}}_5)\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{l}\mathrm{j}\mathrm{n}}=-\mathrm{\rho }{\mathrm{\omega }}^2{\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{k}}\end{array}$$

(41)

$$\begin{array}{l}({\mathrm{c}}_{12}+{\mathrm{c}}_{66})\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{U}}_{\mathrm{l}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{66}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^{\mathrm{\alpha },\mathrm{\epsilon }}}^{\mathrm{x}}{\mathrm{V}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{22}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^2}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{44}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^{\mathrm{\beta }}}^{\mathrm{z}}{\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{n}}\\ +({\mathrm{c}}_{23}+{\mathrm{c}}_{44})\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{m}\mathrm{n}}+({\mathrm{e}}_2+{\mathrm{e}}_4)\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{m}\mathrm{n}}=-\mathrm{\rho }{\mathrm{\omega }}^2{\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{k}}\end{array}$$

(42)

$$\begin{array}{l}({\mathrm{c}}_{13}+{\mathrm{c}}_{55})\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{n}}+({\mathrm{c}}_{23}+{\mathrm{c}}_{44})\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{n}}+({\mathrm{c}}_{55}-{\mathrm{K}}_2)\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^{\mathrm{\alpha },\mathrm{\epsilon }}}^{\mathrm{x}}{\mathrm{W}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+({\mathrm{c}}_{44}-{\mathrm{K}}_2)\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^2}^{\mathrm{y}}{\mathrm{W}}_{\mathrm{i}\mathrm{m}\mathrm{k}}\\ +{\mathrm{c}}_{33}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^{\mathrm{\beta }}}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{e}}_5\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^{\mathrm{\alpha },\mathrm{\epsilon }}}^{\mathrm{x}}{\mathrm{\varphi }}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{e}}_4\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^2}^{\mathrm{y}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{e}}_3\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^{\mathrm{\beta }}}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{K}}_1{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}+{\mathrm{K}}_3{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}^3+{\mathrm{C}}_{\mathrm{d}}(\mathrm{i}\mathrm{\omega }){\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}=-\mathrm{\rho }{\mathrm{\omega }}^2{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}\end{array}$$

(43)

$$\begin{array}{l}({\mathrm{e}}_1+{\mathrm{e}}_5)\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{n}}+({\mathrm{e}}_2+{\mathrm{e}}_4)\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{n}}+{\mathrm{e}}_5\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^{\mathrm{\alpha },\mathrm{\epsilon }}}^{\mathrm{x}}{\mathrm{W}}_{\mathrm{l}\mathrm{j}\mathrm{k}}\\ +{\mathrm{e}}_4\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^2}^{\mathrm{y}}{\mathrm{W}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{e}}_3\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^{\mathrm{\beta }}}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}-{\mathrm{\eta }}_1\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^{\mathrm{\alpha },\mathrm{\epsilon }}}^{\mathrm{x}}{\mathrm{\varphi }}_{\mathrm{l}\mathrm{j}\mathrm{k}}-{\mathrm{\eta }}_2\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^2}^{\mathrm{y}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{m}\mathrm{k}}-{\mathrm{\eta }}_3\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^{\mathrm{\beta }}}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\end{array}$$

(44)

The nonlinear algebraic system may be written as follow:

$$\left[\begin{array}{c}\begin{array}{cc}\begin{array}{cc}\left[{\mathrm{K}}_{11}\right]& \left[{\mathrm{K}}_{12}\right]\end{array}& \begin{array}{cc}\left[{\mathrm{K}}_{13}\right]& \left[{\mathrm{K}}_{14}\right]\end{array}\end{array}\\ \begin{array}{cc}\begin{array}{cc}\begin{array}{c}\left[{\mathrm{K}}_{21}\right]\\ \left[{\mathrm{K}}_{31}\right]\\ \left[{\mathrm{K}}_{41}\right]\end{array}& \begin{array}{c}\left[{\mathrm{K}}_{22}\right]\\ \left[{\mathrm{K}}_{32}\right]\\ \left[{\mathrm{K}}_{42}\right]\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\left[{\mathrm{K}}_{23}\right]\\ \left[{\mathrm{K}}_{33}\right]\\ \left[{\mathrm{K}}_{43}\right]\end{array}& \begin{array}{c}\left[{\mathrm{K}}_{24}\right]\\ \left[{\mathrm{K}}_{34}\right]\\ \left[{\mathrm{K}}_{44}\right]\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{l}\left[\mathrm{U}\right]\\ \left[\mathrm{V}\right]\\ \left[\mathrm{W}\right]\\ \left[\mathrm{\varphi }\right]\end{array}\right]=\mathrm{\omega }\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}-{\mathrm{C}}_{\mathrm{d}}(\mathrm{i})\left[\mathrm{I}\right]\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{l}\left[\mathrm{U}\right]\\ \left[\mathrm{V}\right]\\ \left[\mathrm{W}\right]\\ \left[\mathrm{\varphi }\right]\end{array}\right]+{\mathrm{\omega }}^2\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{c}-\mathrm{\rho }\left[\mathrm{I}\right]\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ -\mathrm{\rho }\left[\mathrm{I}\right]\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}-\mathrm{\rho }\left[\mathrm{I}\right]\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{l}\left[\mathrm{U}\right]\\ \left[\mathrm{V}\right]\\ \left[\mathrm{W}\right]\\ \left[\mathrm{\varphi }\right]\end{array}\right]$$

(45)

 

where, [I] is the unit matrix.

Also, this system can be considered as follow:

$$[\mathrm{K}][\mathrm{X}]=\mathrm{\omega }[\mathrm{R}][\mathrm{X}]+{\mathrm{\omega }}^2[\mathrm{H}][\mathrm{X}]$$

(46)

Then, the overall nonlinear algebraic system for all layers with continuity conditions can be put in the form:

$$\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}{[\mathrm{K}]}^1\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ {[\mathrm{K}]}^2\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}\ddots \\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ {[\mathrm{K}]}^{\mathrm{L}}\end{array}\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{l}{\left[\mathrm{X}\right]}^1\\ {\left[\mathrm{X}\right]}^2\\ ⋮\\ {\left[\mathrm{X}\right]}^{\mathrm{L}}\end{array}\right]=\mathrm{\omega }\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}{[\mathrm{R}]}^1\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ {[\mathrm{R}]}^2\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}\ddots \\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ {[\mathrm{R}]}^{\mathrm{L}}\end{array}\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{l}{\left[\mathrm{X}\right]}^1\\ {\left[\mathrm{X}\right]}^2\\ ⋮\\ {\left[\mathrm{X}\right]}^{\mathrm{L}}\end{array}\right]+{\mathrm{\omega }}^2\left[\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\left[\mathrm{H}\right]}^1\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ {\left[\mathrm{H}\right]}^2\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}\ddots \\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ {\left[\mathrm{H}\right]}^{\mathrm{L}}\end{array}\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{l}{\left[\mathrm{X}\right]}^1\\ {\left[\mathrm{X}\right]}^2\\ ⋮\\ {\left[\mathrm{X}\right]}^{\mathrm{L}}\end{array}\right]$$

(47)

Where, L denotes to the number of layers.

Boundary conditions may also be explained as follow

For the simply supported end:

$$\begin{array}{l}{\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}}=0,\\ {\mathrm{\sigma }}_{\mathrm{x}}=0\left({\mathrm{c}}_{11}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{12}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{13}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{e}}_1\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\right)\mathrm{at}\mathrm{x}=0,\mathrm{a}\end{array}$$

(48)

$$\begin{array}{l}{\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}}=0,\\ {\mathrm{\sigma }}_{\mathrm{y}}=0\left({\mathrm{c}}_{12}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{22}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{23}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{e}}_2\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\right)\mathrm{at}\mathrm{y}=0,\mathrm{b}\end{array}$$

(49)

$$\begin{array}{l}\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{d}\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{d}:{\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}}=0,\mathrm{at}\mathrm{x}=0,\mathrm{a}\\ {\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}}={\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}}=0,\mathrm{a}\mathrm{t}\mathrm{y}=0,\mathrm{b}\end{array}$$

(50)

For the free supported end:

$$\begin{array}{l}{\mathrm{\sigma }}_{\mathrm{x}}=0\left({\mathrm{c}}_{11}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{12}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{13}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{e}}_1\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\right)\\ {\mathrm{\tau }}_{\mathrm{x}\mathrm{z}}=0\left({\mathrm{c}}_{55}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{c}}_{55}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{W}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{e}}_5\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{\varphi }}_{\mathrm{l}\mathrm{j}\mathrm{k}}=0\right)\\ {\mathrm{\tau }}_{\mathrm{x}\mathrm{y}}=0\left({\mathrm{c}}_{66}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{U}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{66}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{V}}_{\mathrm{l}\mathrm{j}\mathrm{k}}=0\right)\\ {\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}}=0,\mathrm{at}\mathrm{x}=0,\mathrm{a}\end{array}$$

(51)

$$\begin{array}{l}{\mathrm{\sigma }}_{\mathrm{y}}=0\left({\mathrm{c}}_{12}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{22}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{23}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{e}}_2\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\right)\\ {\mathrm{\tau }}_{\mathrm{y}\mathrm{z}}=0\left({\mathrm{c}}_{44}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{c}}_{44}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{W}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{e}}_4\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{m}\mathrm{k}}=0\right)\\ {\mathrm{\tau }}_{\mathrm{x}\mathrm{y}}=0\left({\mathrm{c}}_{66}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{U}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{66}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{V}}_{\mathrm{l}\mathrm{j}\mathrm{k}}=0\right)\\ {\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}}=0,\mathrm{at}\mathrm{y}=0,\mathrm{b}\end{array}$$

(52)

For the top and the bottom surfaces, the boundary condition may be written as:

For closed circuit (grounded):

$$\begin{array}{l}{\mathrm{\sigma }}_{\mathrm{z}}=0\left({\mathrm{c}}_{13}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{23}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{33}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{e}}_3\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\right),\\ {\mathrm{\tau }}_{\mathrm{x}\mathrm{z}}=0\left({\mathrm{c}}_{55}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{c}}_{55}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{W}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{e}}_5\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{\varphi }}_{\mathrm{l}\mathrm{j}\mathrm{k}}=0\right),\\ {\mathrm{\tau }}_{\mathrm{y}\mathrm{z}}=0\left({\mathrm{c}}_{44}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{c}}_{44}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{W}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{e}}_4\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{m}\mathrm{k}}=0\right),\\ {\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}}=0\mathrm{at}\mathrm{z}=0,\mathrm{h}\end{array}$$

(53)

For open circuit:

$$\begin{array}{l}{\mathrm{\sigma }}_{\mathrm{z}}=0\left({\mathrm{c}}_{13}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{c}}_{23}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{c}}_{33}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{e}}_3\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\right),\\ {\mathrm{\tau }}_{\mathrm{x}\mathrm{z}}=0\left({\mathrm{c}}_{55}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{c}}_{55}\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{W}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{e}}_5\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{\varphi }}_{\mathrm{l}\mathrm{j}\mathrm{k}}=0\right),\\ {\mathrm{\tau }}_{\mathrm{y}\mathrm{z}}=0\left({\mathrm{c}}_{44}\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{n}}+{\mathrm{c}}_{44}\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{W}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{e}}_4\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{m}\mathrm{k}}=0\right),\\ {\mathrm{D}}_{\mathrm{z}}=0\left({\mathrm{e}}_1\sum _{\mathrm{l}=1}^{{\mathrm{N}}_{\mathrm{x}}}{{\mathrm{C}}_{\mathrm{i}\mathrm{l}}^1}^{\mathrm{x}}{\mathrm{U}}_{\mathrm{l}\mathrm{j}\mathrm{k}}+{\mathrm{e}}_2\sum _{\mathrm{m}=1}^{{\mathrm{N}}_{\mathrm{y}}}{{\mathrm{C}}_{\mathrm{j}\mathrm{m}}^1}^{\mathrm{y}}{\mathrm{V}}_{\mathrm{i}\mathrm{m}\mathrm{k}}+{\mathrm{e}}_3\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{n}}-{\mathrm{\eta }}_3\sum _{\mathrm{n}=1}^{{\mathrm{N}}_{\mathrm{z}}}{{\mathrm{C}}_{\mathrm{k}\mathrm{n}}^1}^{\mathrm{z}}{\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{n}}=0\right)\\ \mathrm{at}\mathrm{z}=0,\mathrm{h}\end{array}$$

(54)

The continuity conditions at the interface locations can be demonstrated as follow:

$$\begin{array}{l}{({\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({\mathrm{U}}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1},{({\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1},{({\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({\mathrm{W}}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1},{({\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({\mathrm{\varphi }}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1},\\ {({({\mathrm{\sigma }}_{\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({({\mathrm{\sigma }}_{\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1},{({({\mathrm{\tau }}_{\mathrm{x}\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({({\mathrm{\tau }}_{\mathrm{x}\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1},{({({\mathrm{\tau }}_{\mathrm{y}\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({({\mathrm{\tau }}_{\mathrm{y}\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1},{({({\mathrm{D}}_{\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}}={({({\mathrm{D}}_{\mathrm{z}})}_{\mathrm{i}\mathrm{j}\mathrm{k}})}_{\mathrm{p}+1}\end{array}$$

(55)

Where p = (1, L−1), L is the number of layers.

Numerical results and discussions

A parametric study is introduced to discuss the influence of a/h ratio, delamination length, delamination position, damper modulus parameter, and boundary conditions on the obtained fundamental frequency.

For practical purpose, the field quantities are put in dimensionless terms as ¹⁰ :

$$\begin{array}{l}(U^{*},V^{*},W^{*})=\frac{(U,V,W)}{h},x^{*}=\frac{x}{a},\\ y^{*}=\frac{y}{b},z^{*}=\frac{z}{h},C_{ij}^{*}=\frac{C_{ij}}{E_2},\\ ({{\sigma }_x}^{*},{{\sigma }_y}^{*},{{\sigma }_z}^{*},{{\tau }_{xy}}^{*},{{\tau }_{xz}}^{*},{{\tau }_{yz}}^{*})=\frac{({\sigma }_x,{\sigma }_y,{\sigma }_z,{\tau }_{xy},{\tau }_{xz},{\tau }_{yz})}{E_2},\\ {\overline{C}}_d=\frac{C_d{a}^2}{\sqrt{\mathrm{\rho }hD}},D=\frac{E_c{h}^3}{12(1-{\upsilon }^2)}\end{array}$$

(56)

where, $(U^{*},V^{*},W^{*})$ represent normalized displacements, $({{\sigma }_x}^{*},{{\sigma }_y}^{*},{{\sigma }_z}^{*})$ refer to the normalized stresses, $({\tau }_{xy},{\tau }_{xz},{\tau }_{yz})$ are normalized shear stresses, $(D)$ refer to the flexural rigidity for the bottom layer of composite, and ${\overline{C}}_d$ is normalized damper modulus parameter.

The dimensionless frequencies $\mathrm{\Omega }\mathrm{and}\overline{\omega }$ are obtained as ¹⁰ :

$$\mathrm{\Omega }=\omega \mathrm{h}\sqrt{\frac{{\mathrm{\rho }}_{\mathrm{c}}}{E_2}},\overline{\omega }=\mathrm{\Omega }{\left(\frac{a}{h}\right)}^2$$

(57)

Numerical quadrature schemes are designed, for each technique of DQM techniques, to solve free vibration problems involving delaminated piezoelectric composite plates that rest on a viscoelastic foundation. A comparative study is presented to compare the obtained results with these schemes against the published exact and numerical solutions to check the validity of these schemes.

Plate (a) is composed of two composite layers with orientation angles $[0^{\mathrm{o}}/{90}^{\mathrm{o}}]$ . There is a delamination at interface between the layers, and this plate rests on a viscoelastic foundation, as shown in Figure 2. The material properties of the layers are $E_1=40E_2,E_3=E_2,G_{12}=G_{13}=0.6E_2,G_{23}=0.5E_2,{\mathrm{\nu }}_{12}={\mathrm{\nu }}_{23}={\mathrm{\nu }}_{13}=0.25$ . Tables 1 and 2 show a convergence study for the four techniques of DQM to determine the suitable technique, grid size $({\mathrm{N}}_{\mathrm{x}}\times {\mathrm{N}}_{\mathrm{y}}\times {\mathrm{N}}_{\mathrm{z}})$ for each layer, bandwidth (2M + 1), and regularization parameter $\sigma$ that give good accurate fundamental frequency parameter $\overline{\omega }$ of this plate with simply supported on all four edges at a/b = 1, a/h = 5, ${\overline{C}}_d$  = 0, ${\mathrm{a}}_1={\mathrm{a}}_3$ , and delamination length ${\mathrm{a}}_2$  = 0. Table 1 presents the results coming from PDQM, Sinc DQM, and DLK DQM with different grid sizes. Also, different fractional order derivatives were applied. From this table, one can conclude that these PDQM gives accurate and convergent results. PDQM and DLK DQM results agree well with Chen and Lü ¹ and Matsunaga ³ at grid size 7 × 5 × 5 but Sinc DQM method is inaccurate at the studied grid size. Also, DLK DQM gives accurate results at grid size 7 × 5 × 5 and M = 1. Further, CPU time is calculated in Table 1. The execution time for DLK DQM is less than PDQM. Moreover, Table 1 presents the effect of fractional derivative α on the obtained results. Fundamental frequency and CPU time at α = 1.98 is more than at α = 2. Table 2 presents results of RSK DQM at different grid sizes, different values of bandwidth parameter (M) with values from 1 up to 5, and regularization parameter $\mathrm{\sigma }$ with values from 0.0335 h to 0.505 h. Also, this table proves that results of RSKDQM at different fractional derivative β. Results are more accurate, efficient and convergent more than other proposed techniques. So, one may use this technique with grid size $7\times 6\times 5$ , bandwidth parameter (M) equal 1, and regularization parameter $\mathrm{\sigma }$ equal 0.0355 h to get accurate result that is near to the previous ones.

Figure 2.

Configurations of laminated plates (a) which is rested on viscoelastic foundation.

Table 1.

Effect of grid size and bandwidth on fundamental frequency $\overline{\omega }$ , for SSSS plate (a) using PDQM, Sinc DQM and DLK DQM (a/b = 1, a/h = 5, ${\overline{C}}_d$  =  $a_2$  = 0 and M = 1).

Grid size	PDQM (α = 2)		PDQM (α = 1.98)		Sinc DQM (α = 2)		Sinc DQM (α = 1.98)		DLK DQM (α = 2)		DLK DQM (α = 1.98)
	$\overline{\omega }$	CPU time (sec)	$\overline{\omega }$	CPU time (sec)	$\overline{\omega }$	CPU time (sec)	$\overline{\omega }$	CPU time (sec)	$\overline{\omega }$	CPU time (sec)	$\overline{\omega }$	CPU time (sec)
5 × 5 × 5	8.3940	1.35	12.1392	2.182	84.2472	0.971	85.2537	1.32	9.4868	0.90	9.4868	1.33
6 × 5 × 5	8.4325	1.36	8.8566	2.869	32.1889	1.49	32.8662	1.94	9.4868	1.52	9.4868	2.00
7 × 5 × 5	8.4391	2.78	8.9063	3.5131	56.1408	2.46	56.0680	3.16	8.8568	2.45	9.4868	3.25
8 × 5 × 5	8.4383	3.94	8.7299	4.898	35.9359	4.02	40.8973	4.67	9.4868	4.13	9.4868	4.96
9 × 5 × 5	8.4383	6.02	8.7608	7.608	43.8864	6.22	49.4084	6.98	9.3554	6.20	9.4868	7.05
Chen and Lü 1	8.5269
Matsunaga 3	8.5301

Table 2.

Effect of bandwidth, and regularization parameter $\sigma$ on obtained fundamental frequency $\overline{\omega }$ , for SSSS plate (a) using RSK DQM with grid size $7\times 6\times 5$ (a/b = 1, a/h = 5, ${\overline{C}}_d$  =  $a_2$  = 0).

β	$\sigma$	M = 1	M = 2	M = 3	M = 4	M = 5
	0.0335 h	8.5225	8.5225	8.5225	8.5225	8.5225
	0.0355 h	8.5261	8.5261	8.5261	8.5261	8.5261
2	0.0478 h	8.5261	8.5261	8.5261	8.5261	8.5261
	0.0481 h	8.5261	8.5261	8.5261	8.5261	8.5261
	0.0499 h	8.5261	8.5261	8.5261	8.5261	8.5261
	0.0505 h	8.5261	8.5261	8.5261	8.5261	8.5261
	Execution time	3.696958 sec at M = 1, and $\sigma$  = 0.0355 h
	0.0335 h	9.9358	9.9358	9.9358	9.9358	9.9358
	0.0355 h	9.6342	9.6342	9.6342	9.6342	9.6342
	0.0478 h	9.6342	9.6342	9.6342	9.6342	9.6342
1.98	0.0481 h	9.6342	9.6342	9.6342	9.6342	9.6342
	0.0499 h	9.6342	9.6342	9.6342	9.6342	9.6342
	0.0505 h	9.6342	9.6342	9.6342	9.6342	9.6342
	Execution time	5.708946 sec at M = 1, and $\sigma$  = 0.0355 h
Chen and Lü 1	8.5269
Matsunaga 3	8.5301

It is noted from previous tables that, RSK DQM give the best accurate and efficient results compared to the other techniques, So RSK DQM method is employed for deducing the parametric studies. Tables 3 and 4 show the obtained fundamental frequency parameter $\overline{\mathrm{\omega }}$ of this plate with ratio a/b = 1 under different boundary conditions, different values of a/h ratio, and different delamination length ratios ${\mathrm{a}}_2/\mathrm{a}$ . The obtained results have a good convergent and are similar to the compared ones. Table 3 displays fundamental frequency parameter $\overline{\mathrm{\omega }}$ of the plate (a) at different fractional derivatives ε with a delamination at middle of the x-direction span. The obtained results of SSSS plate with zero delamination length are compared with the results of Chen and Lü. ¹ It is concluded from previous tables, that fundamental frequency increases with decreasing fractional derivatives α and β. Also, fundamental frequency decreases with decreasing fractional derivatives ε. In addition, Table 4 illustrates the fundamental frequency parameter of the plate (a) $[0^{°}/{90}^{°}]$ with a delamination where its center is located at the middle of the x-direction span $({\mathrm{a}}_1={\mathrm{a}}_3)$ . In addition, the used values of damper modulus parameter ${\overline{C}}_d$ are 0, 0.5, 1, 1.5, 2, and 2.5. RSKDQM is used with a grid size of ${\mathrm{N}}_{\mathrm{x}}\times {\mathrm{N}}_{\mathrm{y}}\times {\mathrm{N}}_{\mathrm{z}}=(7\times 6\times 5)$ for each layer. The assumed values of length ratio of the delamination ${\mathrm{a}}_2/\mathrm{a}$ are 0, 0.3, 0.6, and 0.9. The obtained results for the fundamental frequency parameter of SSSS plate, simply supported on all four edges, are appeared in Table 4.

Table 3.

Comparison of the fundamental frequency $\overline{\omega }$ , for SSSS plate (a): $[0^o/{90}^o]$ that has a central delamination using RSK DQM at M = 1, and $\sigma$  = 0.0355 h with a grid size $7\times 6\times 5$ (a/b = 1, ${\overline{C}}_d$  = 0, $a_1=a_3$ ).

ε	a/h	Without delamination 1	Delamination length ratio $a_2/a$
			0.05	0.1	0.25	0.5	0.75	0.9
1	5	8.5225	8.3450	8.3450	8.2405	7.6647	7.0847	6.5670
	10	10.33645	10.2981	10.2981	10.0615	9.7318	9.4741	9.3576
	20	11.03678	10.8295	10.8295	10.8295	10.4751	9.8020	9.8020
	25	11.13208	10.9194	10.9194	10.8032	10.4092	10.0005	9.7387
	50	11.26359	11.0396	11.0396	11.0396	10.6841	9.9423	9.3120
	100	11.29729	11.0693	11.0693	11.0693	10.7142	9.9599	9.9599
0.99	5	8.3881	8.2537	8.2537	8.2537	8.0022	7.6791	7.1940
	10	10.0147	9.8568	9.8568	9.8568	9.5174	8.9477	8.3856
	20	9.9286	9.7535	9.7535	9.7535	9.3621	8.3182	7.2983
	25	9.3987	9.2120	9.2120	9.2120	8.7928	7.1432	3.1102
	50	3.2840	3.7744	3.7744	3.7744	4.6649	5.4985	3.8524
	100	14.9233	17.5713	17.5713	17.5713	8.1152	7.5451	20.3600

Table 4.

The fundamental frequency of SSSS plate (a): [0/90] with central delamination crack at the interface of the two layers. a/b = 1, a/h = 5, N = 7 × 6 × 5 using RSKDQM. E1 = 40E2, E3 = E2, G12 = G13 = 0.6E2, G23 = 0.5E2, v12 = v13 = v23 = 0.25, a/b = 1 (a1 = a3).

Cd	Delamination crack width a2/a
	0	0.3	0.6	0.9
0	8.5295	8.3999	7.8486	7.3683
0.5	8.5292 ± 0.0710i	8.3996 ± 0.0717i	7.8482 ± 0.0735i	7.3679 ± 0.0744i
1	8.5283 ± 0.1420i	8.3987 ± 0.1434i	7.8472 ± 0.1470i	7.3668 ± 0.1488i
1.5	8.5269 ± 0.2130i	8.3972 ± 0.2152i	7.8455 ± 0.2206i	7.3649 ± 0.2232i
2	8.5248 ± 0.2840i	8.3951 ± 0.2869i	7.8431 ± 0.2941i	7.3623 ± 0.2975i
2.5	8.5222 ± 0.3550i	8.3923 ± 0.3586i	7.8400 ± 0.3676i	7.3589 ± 0.3719i
Chen and Lü 1 at a2/a = 0	8.5269
Matsunaga 3 at a2/a = 0	8.5301

Furthermore, a parametric study is introduced in Figures 3 and 4 to show the influence of a/h ratio, the normalized damper modulus parameter ${\overline{C}}_d$ , the delamination length ratio ${\mathrm{a}}_2/\mathrm{a}$ , and different fractional derivatives on the frequencies of plate (a). From Figures 3–4, one can conclude that, the natural frequency increases with the increasing of length to thickness ratio a/h for SSSS plates until (a/h) > 15, frequency reached to steady state. In addition, the natural frequency decreases with the increasing of delamination length ratio ${\mathrm{a}}_2/\mathrm{a}$ until $({\mathrm{a}}_2/\mathrm{a})\ge 2$ , frequency results are convergent. Also, the natural frequency decreases with the increasing of the normalized damper modulus parameter ${\overline{C}}_d$ for SSSS plates. Further, fundamental frequency increases with decreasing fractional derivatives α as shown in Figure 3. In addition, fundamental frequency decreases with decreasing fractional derivatives β as shown in Figure 4.

Figure 3.

Variation of obtained fundamental frequency parameter $\overline{\omega }$ with the ratio a/h, for the plate (a) with SSSS boundary conditions using RSK DQM with a grid size $7\times 6\times 5$ , (2M + 1) = 3, ${\overline{C}}_d$  = 0, and $\sigma$  = 0.0355 h for each layer $a_2/a$  = 0, a/b = 1, a/h = 5, $a_1=a_3$ : (a) α = 2, (b) α = 1.99, and (c) α = 1.98.

Figure 4.

Variation of obtained fundamental frequency parameter $\overline{\omega }$ with the ratio a/h, for the plate (a) with SSSS boundary conditions using RSK DQM with a grid size $7\times 6\times 5$ , (2M + 1) = 3, and $\sigma$  = 0.0355 h for each layer $a_2/a$  = 0, a/b = 1, a/h = 5, $a_1=a_3$ : (a) β = 2, ${\overline{C}}_d$  = 1, (b) β = 1.99, ${\overline{C}}_d$  = 1, and (c) β = 1.99, ${\overline{C}}_d$  = 0.

Conclusion

The free vibration of a delaminated fractional piezoelectric composite plate resting on a viscoelastic foundation with SSSS boundary conditions is effectively investigated using a numerical scheme based on DQM. The three-dimensional theory of piezoelectricity and elasticity is used to formulate the governing equations. Furthermore, parametric analysis is used to examine how the resulting fundamental frequency is affected by the a/h ratio, delamination length, delamination position, foundation parameters, and damper modulus parameter. PDQM, RSK DQM, Sinc DQM and DLK DQM techniques are introduced. RSK DQM method at M = 1, and $\mathrm{\sigma }$  = 0.0355 h gives more accurate, convergent and efficient results compared to other proposed techniques. Also, a parametric study is employed using RSK DQM. Form parametric study, one can conclude the following notes:

• The fundamental frequency $\overline{\omega }$ increases with an increasing ratio (a/h) for SSSS plates.

• The fundamental frequency $\overline{\omega }$ decreases with increasing in delamination length ratio ${\mathrm{a}}_2/\mathrm{a}$ for SSSS plates.

• The fundamental frequency $\overline{\omega }$ increases with decreasing fractional derivatives α and β. Also, fundamental frequency decreases with decreasing fractional derivatives ε.