# Vibration Analysis of Piezoelectric Composite Plate Resting on Nonlinear Elastic Foundations Using Sinc and Discrete Singular Convolution Differential Quadrature Techniques

Authors: Ola Ragb; Mohamed Salah; M. S. Matbuly; R. B. M. Amer
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this work, free vibration of the piezoelectric composite plate resting on nonlinear elastic foundations is examined. The three-dimensionality of elasticity theory and piezoelectricity is used to derive the governing equation of motion. By implementing two differential quadrature schemes and applying different boundary conditions, the problem is converted to a nonlinear eigenvalue problem. The perturbation method and iterative quadrature formula are used to solve the obtained equation. Numerical analysis of the proposed schemes is introduced to demonstrate the accuracy and efficiency of the obtained results. The obtained results are compared with available results in the literature, showing excellent agreement. Additionally, the proposed schemes have higher efficiency than previous schemes. Furthermore, a parametric study is introduced to investigate the effect of elastic foundation parameters, different materials of sensors and actuators, and elastic and geometric characteristics of the composite plate on the natural frequencies and mode shapes.


---

## Body

## 1. Introduction

The increase in the use of piezoelectric composite materials resting on elastic foundation structures, especially in aerospace, automotive, and marine environments, leads to the rise of difficulties in nonlinear vibrations in various modern engineering challenges. The elastic foundations play an important role that preserve the structural system under oscillations and avoid mechanical failures. The main role of nonlinear foundations is developing the accuracy of the model to describe the behavior of our system. The cubic nonlinearity is the most common type used in recent research [1–8].Due to the mathematical complexity of such problems, only limited cases can analytically be solved. A number of exact methods of the vibration problems on elastic foundations are issued in [9–14]. Literature on the numerical solution of the research subject is sparse. Typical useful numerical methods, such as the finite element [15–17], the boundary element [18], differential quadrature [19], and state space [20, 21] techniques, are used to solve such problems. One of the disadvantages of these methods is consuming too much time and effort to get the solution [22–26]. Additionally, there is no research in piezoelectric composite plate resting on nonlinear elastic foundations, whether analytically or numerically.The sinc differential quadrature method (SDQM) [27–32] and the discrete singular convolution differential quadrature method (DSCDQM) [33–42] are more reliable than the polynomial-based DQM. The nonlinear partial differential equations are reduced to linear equations using perturbation [43] and iterative quadrature methods [44]. Previous DQ techniques are then employed to reduce the problem to the eigenvalue or bending problem. The same schemes are used for free vibration analysis of piezoelectric nanobeams [45, 46]. Recent studies have investigated the vibration analysis of composite plate resting on elastic foundations, using DSCDQM [47–58].As far as the authors are aware, SDQM and DSCDQM have not been examined for vibration analysis of composite piezoelectric plate materials resting on nonlinear elastic foundations. Based on these versions, numerical schemes are designed for free vibration of piezoelectric composites. The natural frequencies are obtained and compared with previous analytical and numerical frequencies. For each scheme, the convergence and efficiency are verified. Additionally, a parametric study is introduced to investigate the influence of elastic foundation parameters and elastic and geometric characteristics of the composite on the vibrated results.


## 2. Formulation of the Problem

Consider a three-dimensional piezoelectric composite with0≤x≤a,0≤y≤b,0≤z≤h, where a, b, and h are the length, width, and total thickness of the composite, respectively. This composite is polarized in the z direction and consists of m layers with different types of materials. The plate rests on a three-parameter foundation model, as shown in Figure 1.Figure 1
Piezoelectric composite (sensor is PZT-4 and actuator is Ba2NaNb5O15) resting on nonlinear elastic foundations.Based on the theory of elasticity and piezoelectricity, the equations of motion and the charge equation of electrostatic can be written as follows [59, 65–69]:(1)σx,x+τxy,y+τxz,z=ρu,tt,(2)τxy,x+σy,y+τyz,z=ρv,tt,(3)τxz,x+τyz,y+σz,z+k1w+k3w3−k2w,xx+w,yy=ρw,tt,(4)Dx,x+Dy,y+Dz,z=0,where σx,σy,σz, u,v,w, and Dx,Dy,Dz are the stresses, displacement, and induction field in the x, y, and z directions, respectively; τyz,τxz,τxy are the shear stresses; ρ is the density of the material; and (k1, k2, k3) are linear Winkler foundation, linear Pasternak foundation, and nonlinear Winkler foundation parameters, respectively.The relation between mechanical and electric material properties is constitutive equations, which can be written as follows:(5)σxσyσzτyzτxzτxy=C11C12C13000C12C11C13000C13C13C33000000C44000000C55000000C66u,xv,yw,zw,y+v,zw,x+u,zv,x+u,y+0000e50000e40000e3000Tϕ,xϕ,yϕ,z,(6)DxDyDz=0000e50000e40000e3000u,xv,yw,zw,y+v,zw,x+u,zv,x+u,y−η1000η2000η3ϕ,xϕ,yϕ,z,where C, e, and η are the components of the effective elastic, piezoelectric, and dielectric constants of the same piezoelectric material, respectively, and ϕ is the electric potential.For harmonic behavior, one can assume that(7)ux,z,t=Ueiωt,vx,z,t=Veiωt,wx,z,t=Weiωt,ϕx,z,t=Φeiωt,where ω is the natural frequency of the plate and i=−1. U,W,andΦ are the amplitudes for u,w, and ϕ, respectively.The elastic material constants can be determined as follows, using the reciprocal theorem [41]:(8)C11=E1E2−v232E3E2−v122E22/E1−v232E3−v13E2E3/E1v13−2v12v23,(9)C12=C21=v12E2+v13v23E31−v122E2/E1−v232E3/E2−v13E3/E1v13−2v12v23,(10)C13=C31=E3v13+v12v231−v122E2/E1−v232E3/E2−v13E3/E1v13−2v12v23,(11)C22=E2E1−v132E3E1−v122E2−v232E1E3/E2−v13E3v13−2v12v23,(12)C23=C32=E3v23E1+v12v13E2E1−v122E2−v232E1E3/E2−v13E3v13−2v12v23,(13)C33=E3E1−v122E2E1−v122E2−v232E1E3/E2−v13E3v13−2v12v23,(14)C44=G23,C55=G13,C66=G12,where Ep, Gpq, and vpq (p, q = 1, 2, 3) are Young’s moduli, shear moduli, and Poisson’s ratios.The isotropic material constants can be expressed as follows [56]:(15)C11=C12=C13=C22=C23=C31=C32=C33=νΕ1+ν1−2ν+2Ε21+ν,C44=C55=C66=Ε21+ν.Substituting equations (5)–(16) into (1)–(4), the problem can be reduced to a quasistatic one, as follows:(16)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W,xz+e1+e5Φ,xz=−ρω2U,(17)C66+C12U,xy+C662V,xx+C11V,yy+C44V,zz+C13+C44W,yz+e2+e4Φ,yz=−ρω2V,(18)C55+C13U,xz+C13+C44V,yz+C55−k2W,xx+C44−k2W,yy+C33W,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1W+k3W3=−ρω2W,(19)e1+e5U,xz+e2+e4V,yz+e5W,xx+e4W,yy+e3W,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.The boundary conditions can be described as follows:(1)
For a simply supported edge (S):(20)w=v=σx=0,atx=0,a,w=u=σy=0,aty=0,b.(2)
For a clamped edge (C):(21)w=v=u=0,atx=0,ay=0,b.(3)
For a free edge (F):(22)σx=τxz=τxy=0,atx=0,a,σy=τyz=τxy=0,aty=0,b.Mechanical and electrical boundary conditions at the lower and upper surfaces of the composite are as follows:(23)σz=τxz=τzy=Dz=0,atz=0,σz=τxz=τzy=Φ=0,atz=h,To ensure the continuity between electric and elastic layers, the following conditions can be considered:(24)Ux,y,hs−=Ux,y,hs+,Vx,y,hs−=Vx,y,hs+Wx,y,hs−=Wx,y,hs+,Φx,y,hs−=Φx,y,hs+,(25)Ux,y,hs+hc−=Ux,y,hs+hc+,Vx,y,hs+hc−=Vx,y,hs+hc+,Wx,y,hs+hc−=Wx,y,hs+hc+,Φx,y,hs+hc−=Φx,y,hs+hc+.Additionally, the continuity conditions between different elastic materials are as follows:(26)Ux,y,hk−=Ux,y,hk+,Vx,y,hk−=Vx,y,hk+,Wx,y,hk−=Wx,y,hk+,Φx,y,hk−=Φx,y,hk+.


## 3. Method of Solution

Two differential quadrature techniques are employed to reduce the governing equations into a nonlinear eigenvalue problem, as follows.


### 3.1. Sinc Differential Quadrature Method (SDQM)

A cardinal sine function is used as a shape function, such that the unknownψ and its derivatives can be approximated as a weighted linear sum of nodal values, ψi, (i = −N, N), as follows [27–32]:(27)ψxi=∑j=−NNsinπxi−xj/hxπxi−xj/hxψxj,i=−N,N,hx>0,(28)∂ψ∂xx=xi=∑j=−NNAijxψxj,i=−N,N,(29)∂2ψ∂x2x=xi=∑j=−NNBijxψxj,i=−N,N,where ψ denotes U,V,W, and Φ; N is the number of grid points; and hx is the grid size. The weighting coefficients, Aijx and Bijx, can be determined by differentiating (28) as follows:(30)Aijx=−1i−jhxi−j,i≠j,0,i=j,Bijx=−2−1i−jhx2i−j2i≠j,−π23hx2i=j.


### 3.2. Discrete Singular Convolution Differential Quadrature Method (DSCDQM)

A singular convolution can be defined as follows [33–42]:(31)Ft=T∗ηt=∫−∞∞Tt−xηxdx,where Tt−x is a singular kernel.The DSC algorithm can be applied using many types of kernels. These kernels are applied as shape functions, such that the unknownψ and its derivatives are approximated as a weighted linear sum of ψi, (i = −N, N) over a narrow bandwidth (x−xM,x+xM) [33–37].Two kernels of DSC will be employed as follows:(a)
A Delta Lagrange kernel (DLK) can be used as a shape function, such that the unknownψ and its derivatives can be approximated as a weighted linear sum of nodal values, ψi, (i = −N, N), as follows:(32)ψxi=∑j=−MM∏k=−MMxi−xkxi−xj∏j=−M,j≠kMxj−xkψxj,i=−N,N,M≥1,(33)∂ψ∂xx=xi=∑j=−MMAijxψxj,∂2ψ∂x2x=xi=∑j=−MMBijxψxj,i=−N,N,where 2M + 1 is the effective computational bandwidth. AijxandBijx are defined as follows:(34)Aijx=1xi−xj∏k=−M,k≠i,jMxi−xkxj−xk,i≠j,−∑j=−M,j≠i−MAijx,i=j,Bijx=2Aijx.Aiix−Aijxxi−xj,i≠j,−∑j=−M,j≠iMBijx,i=j.(b)
A regularized Shannon kernel (RSK) can also be used as a shape function, such that the unknownψ and its derivatives can be approximated as a weighted linear sum of nodal values, ψi, (i = −N, N), as follows:(35)ψxi=∑j=−MMsinπxi−xj/hxπxi−xj/hxe−xi−xj2/2σ2ψxj,i=−N,N,σ=r∗hx>0,(36)∂ψ∂xx=xi=∑j=−MMAijxψxj,∂2ψ∂x2x=xi=∑j=−MMBijxψxj,i=−N,N,where σ is a regularization parameter and r is a computational parameter. The weighting coefficients Aijx and Bijx can be defined as follows [33–42]:(37)Aijx=−1i−jhxi−je−hxi−j2/2σ2,i≠j,0,i=j,Bijx=2−1i−j+1hx2i−j2+1σ2e−hx2i−j2/2σ2,i≠j,−1σ2−π23hx2,i=j.Similarly, one can approximateψy,ψz, ψyy,andψzz and calculate Aijy,Aijz, Bijy,andBijz. Two methods are used to transform the governing equation to linear equations.


#### 3.2.1. Perturbation Method

Assuming the solution of equations (17)–(27) to be a power series that can be solved by a perturbation method of the second order [43] by substituting W=Wo+εW1+ε2W2+⋯ into (17)–(27) and equating the terms with the identical powers of ε, where ε is a perturbation parameter:(38)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55Wo,xz+εW1,xz+ε2W2,xz+e1+e5Φ,xz=−ρω2U,(39)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44Wo,yz+εW1,yz+ε2W2,yz+e2+e4Φ,yz=−ρω2V,(40)C55+C13U,xz+C13+C44V,yz+C55−k2Wo,xx+εW1,xx+ε2W2,xx+C44−k2Wo,yy+εW1,yy+ε2W2,yy+C33Wo,zz+εW1,zz+ε2W2,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wo+εW1+ε2W2+k3εWo+εW1+ε2W23=−ρω2Wo+εW1+ε2W2,(41)e1+e5U,xz+e2+e4V,yz+e5Wo,xx+εW1,xx+ε2W2,xx+e4Wo,yy+εW1,yy+ε2W2,yy+e3Wo,zz+εW1,zz+ε2W2,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.Equate the terms ofε0in (39)–(42):(42)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W0,xz+e1+e5Φ,xz=−ρω2U,(43)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W0,yz+e2+e4Φ,yz=−ρω2V,(44)C55+C13U,xz+C13+C44V,yz+C55−k2Wo,xx+C44−k2Wo,yy+C33Wo,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wo=−ρω2Wo,(45)e1+e5U,xz+e2+e4V,yz+e5W0,xx+e4W0,yy+e3W0,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.The previous system is an eigenvalue problem, which is solved to obtain the natural frequencies and the bending deflectionW0.Equate the terms ofε1 in (39)–(42):(46)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W1,xz+e1+e5Φ,xz=−ρω2U,(47)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W1,yz+e2+e4Φ,yz=−ρω2V,(48)C55+C13U,xz+C13+C44V,yz+C55−k2W1,xx+C44−k2W1,yy+C33W1,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1+ρω2W1=−k3Wo3,(49)e1+e5U,xz+e2+e4V,yz+e5W1,xx+e4W1,yy+e3W1,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0,The previous system is a bending problem, which is solved to obtainW1.Equate the terms ofε2 in (39)–(42):(50)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W2,xz+e1+e5Φ,xz=−ρω2U,(51)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W2,yz+e2+e4Φ,yz=−ρω2V,(52)C55+C13U,xz+C13+C44V,yz+C55−k2W2,xx+C44−k2W2,yy+C33W2,zz+e5Φ,xx+e4Φ,yy+e3Φ,zzk1+ρω2W2=−3k3Wo2W1,(53)e1+e5U,xz+e2+e4V,yz+e5W2,xx+e4W2,yy+e3W2,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.The previous system is a bending problem, which is solved to obtainW2.Finally, the series solution can be written as follows:(54)Wnumerical=limε⟶1W=Wo+W1+W2.The convergence condition [43] due to the perturbation method is set as follows:(55)Wi+1Wi<1,wherei=0,1.


#### 3.2.2. Iterative Quadrature

We solved the following iterative system:(56)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55Wi+1,xz+e1+e5Φ,xz=−ρω2U,i=0,1,2,(57)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44Wi+1,yz+e2+e4Φ,yz=−ρω2V,i=0,1,2,(58)C55+C13U,xz+C13+C44V,yz+C55−k2Wi+1,xx+C44−k2Wi+1,yy+C33Wi+1,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wi+1+k3Wi2Wi+1=−ρω2Wi+1,i=0,1,2,(59)e1+e5U,xz+e2+e4V,yz+e5Wi+1,xx+e4Wi+1,yy+e3Wi+1,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0,i=0,1,2.To getWo, we solved the following eigenvalue problem:(60)C11∑l=−NxNxBilxUljk+C66∑m=−NyNyBjmyUimk+C55∑n=−NzNzBknzUijn+C12+C66∑l=−NxNxAilx∑m=−NyNyAjmyVlmk+C13+C55∑l=−NxNxAilx∑n=−NzNzAknzW0,ljn+e1+e5∑l=−NxNxAilx∑n=−NzNzAknzΦljn=−ρω2U,(61)C12+C66∑l=−NxNxAilx∑m=−NyNyAjmyUlmk+C66∑l=−NxNxBilxVljk+C11∑m=−NyNyBjmyVimk+C44∑n=−NzNzBknzVijn+C13+C44∑m=−NyNyAjmy∑n=−NzNzAknzW0,imn+e2+e4∑m=−NyNyAjmy∑n=−NzNzAknzΦimn=−ρω2V,(62)C55+C13∑l=−NxNxAilx∑n=−NzNzAknzUljn+C13+C44∑m=−NyNyAjmy∑n=−NzNzAknzVimn+C55−k2∑l=−NxNxBilxW0,ljk+C44−k2∑m=−NyNyBjmyW0,imk+C33∑n=−NzNzBknzWijn+e5∑l=−NxNxBilxΦljk+e4∑m=−NyNyBjmyΦimk+e3∑n=−NzNzBknzΦijn+k1W0,ijk=−ρω2W0,(63)e1+e5∑l=−NxNxAilx∑n=−NzNzAknzUljn+e2+e4∑m=−NyNyAjmy∑n=−NzNzAknzVimn+e5∑l=−NxNxBilxW0,ljk+e4∑m=−NyNyBjmyW0,imk+e3∑n=−NzNzBknzW0,ijn−η1∑l=−NxNxBilxΦljk−η2∑m=−NyNyBjmyΦimk−η3∑n=−NzNzBknzΦljn=0.The boundary conditions (21)–(27) can also be approximated using two DQMs as follows:(1)
Simply supported (S):(64)Wijk=Vijk=C11∑l=−NxNxAilxUljk+C12∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e1∑n=−NzNzAknzΦijn=0,atx=0,a,Wijk=Uijk=C12∑l=−NxNxAilxUljk+C11∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e2∑n=−NzNzAknzΦijn=0,aty=0,b.(2)
Clamped (C):(65)Wijk=Vijk=Uijk=0,atx=0,ay=0,b.(3)
Free surface (F):(66)C11∑l=−NxNxAilxUljk+C12∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e1∑n=−NzNzAknzΦijn=0,C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=−NxNxAilxΦljk=C66∑l=−NxNxAilxVljk+∑m=−NyNyAjmyUimk=0,atx=0,a,C12∑l=−NxNxAilxUljk+C11∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e2∑n=−NzNzAknzΦijn=0,C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+e4∑m=−NyNyAjmyΦimk=C66∑l=−NxNxAilxVljk+∑m=−NyNyAjmyUimk=0,aty=0,b,Mechanical and electrical boundary conditions at the lower and upper surfaces of the composite are as follows:(67)C13∑l=−NxNxAilxUljk+C13∑m=−NyNyAjmyVimk+C33∑n=−NzNzAknzWijn+e3∑n=−NzNzAknzΦijn=C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=−NxNxAilxΦljk=0C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+e4∑m=−NyNyAjmyΦimk=e1∑l=−NxNxAilxUljk+e2∑m=−NyNyAjmyVimk+e3∑n=−NzNzAknzWijn−η3∑n=−NzNzAknzΦijn=0,atz=0,C13∑l=−NxNxAilxUljk+C13∑m=−NyNyAjmyVimk+C33∑n=−NzNzAknzWijn+e3∑n=−NzNzAknzΦijn=C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=1NxAilxΦljk=0C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+∑m=−NyNyAjmyΦimk=Φijk=0,atz=h.The continuity conditions between the interfaces of layers can be assumed as follows:(68)Wijk1=Wijk2,Vijk1=Vijk2,Uijk1=Uijk2,Φijk1=Φijk2.We have solved the generalized eigenvalue problem [65]:(69)KX=ω2MX,where K is the coefficient matrix of the previous system, M is the mass matrix and can be diagonal with 0 diagonal elements, and ω is the free vibration frequencies squared.For nontrivial solutions for (70), the following determinant should be 0:(70)K−ω2M=0.Equation (71) gives the natural frequencies for the composite plate.


## 3.1. Sinc Differential Quadrature Method (SDQM)

A cardinal sine function is used as a shape function, such that the unknownψ and its derivatives can be approximated as a weighted linear sum of nodal values, ψi, (i = −N, N), as follows [27–32]:(27)ψxi=∑j=−NNsinπxi−xj/hxπxi−xj/hxψxj,i=−N,N,hx>0,(28)∂ψ∂xx=xi=∑j=−NNAijxψxj,i=−N,N,(29)∂2ψ∂x2x=xi=∑j=−NNBijxψxj,i=−N,N,where ψ denotes U,V,W, and Φ; N is the number of grid points; and hx is the grid size. The weighting coefficients, Aijx and Bijx, can be determined by differentiating (28) as follows:(30)Aijx=−1i−jhxi−j,i≠j,0,i=j,Bijx=−2−1i−jhx2i−j2i≠j,−π23hx2i=j.


## 3.2. Discrete Singular Convolution Differential Quadrature Method (DSCDQM)

A singular convolution can be defined as follows [33–42]:(31)Ft=T∗ηt=∫−∞∞Tt−xηxdx,where Tt−x is a singular kernel.The DSC algorithm can be applied using many types of kernels. These kernels are applied as shape functions, such that the unknownψ and its derivatives are approximated as a weighted linear sum of ψi, (i = −N, N) over a narrow bandwidth (x−xM,x+xM) [33–37].Two kernels of DSC will be employed as follows:(a)
A Delta Lagrange kernel (DLK) can be used as a shape function, such that the unknownψ and its derivatives can be approximated as a weighted linear sum of nodal values, ψi, (i = −N, N), as follows:(32)ψxi=∑j=−MM∏k=−MMxi−xkxi−xj∏j=−M,j≠kMxj−xkψxj,i=−N,N,M≥1,(33)∂ψ∂xx=xi=∑j=−MMAijxψxj,∂2ψ∂x2x=xi=∑j=−MMBijxψxj,i=−N,N,where 2M + 1 is the effective computational bandwidth. AijxandBijx are defined as follows:(34)Aijx=1xi−xj∏k=−M,k≠i,jMxi−xkxj−xk,i≠j,−∑j=−M,j≠i−MAijx,i=j,Bijx=2Aijx.Aiix−Aijxxi−xj,i≠j,−∑j=−M,j≠iMBijx,i=j.(b)
A regularized Shannon kernel (RSK) can also be used as a shape function, such that the unknownψ and its derivatives can be approximated as a weighted linear sum of nodal values, ψi, (i = −N, N), as follows:(35)ψxi=∑j=−MMsinπxi−xj/hxπxi−xj/hxe−xi−xj2/2σ2ψxj,i=−N,N,σ=r∗hx>0,(36)∂ψ∂xx=xi=∑j=−MMAijxψxj,∂2ψ∂x2x=xi=∑j=−MMBijxψxj,i=−N,N,where σ is a regularization parameter and r is a computational parameter. The weighting coefficients Aijx and Bijx can be defined as follows [33–42]:(37)Aijx=−1i−jhxi−je−hxi−j2/2σ2,i≠j,0,i=j,Bijx=2−1i−j+1hx2i−j2+1σ2e−hx2i−j2/2σ2,i≠j,−1σ2−π23hx2,i=j.Similarly, one can approximateψy,ψz, ψyy,andψzz and calculate Aijy,Aijz, Bijy,andBijz. Two methods are used to transform the governing equation to linear equations.


### 3.2.1. Perturbation Method

Assuming the solution of equations (17)–(27) to be a power series that can be solved by a perturbation method of the second order [43] by substituting W=Wo+εW1+ε2W2+⋯ into (17)–(27) and equating the terms with the identical powers of ε, where ε is a perturbation parameter:(38)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55Wo,xz+εW1,xz+ε2W2,xz+e1+e5Φ,xz=−ρω2U,(39)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44Wo,yz+εW1,yz+ε2W2,yz+e2+e4Φ,yz=−ρω2V,(40)C55+C13U,xz+C13+C44V,yz+C55−k2Wo,xx+εW1,xx+ε2W2,xx+C44−k2Wo,yy+εW1,yy+ε2W2,yy+C33Wo,zz+εW1,zz+ε2W2,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wo+εW1+ε2W2+k3εWo+εW1+ε2W23=−ρω2Wo+εW1+ε2W2,(41)e1+e5U,xz+e2+e4V,yz+e5Wo,xx+εW1,xx+ε2W2,xx+e4Wo,yy+εW1,yy+ε2W2,yy+e3Wo,zz+εW1,zz+ε2W2,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.Equate the terms ofε0in (39)–(42):(42)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W0,xz+e1+e5Φ,xz=−ρω2U,(43)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W0,yz+e2+e4Φ,yz=−ρω2V,(44)C55+C13U,xz+C13+C44V,yz+C55−k2Wo,xx+C44−k2Wo,yy+C33Wo,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wo=−ρω2Wo,(45)e1+e5U,xz+e2+e4V,yz+e5W0,xx+e4W0,yy+e3W0,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.The previous system is an eigenvalue problem, which is solved to obtain the natural frequencies and the bending deflectionW0.Equate the terms ofε1 in (39)–(42):(46)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W1,xz+e1+e5Φ,xz=−ρω2U,(47)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W1,yz+e2+e4Φ,yz=−ρω2V,(48)C55+C13U,xz+C13+C44V,yz+C55−k2W1,xx+C44−k2W1,yy+C33W1,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1+ρω2W1=−k3Wo3,(49)e1+e5U,xz+e2+e4V,yz+e5W1,xx+e4W1,yy+e3W1,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0,The previous system is a bending problem, which is solved to obtainW1.Equate the terms ofε2 in (39)–(42):(50)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W2,xz+e1+e5Φ,xz=−ρω2U,(51)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W2,yz+e2+e4Φ,yz=−ρω2V,(52)C55+C13U,xz+C13+C44V,yz+C55−k2W2,xx+C44−k2W2,yy+C33W2,zz+e5Φ,xx+e4Φ,yy+e3Φ,zzk1+ρω2W2=−3k3Wo2W1,(53)e1+e5U,xz+e2+e4V,yz+e5W2,xx+e4W2,yy+e3W2,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.The previous system is a bending problem, which is solved to obtainW2.Finally, the series solution can be written as follows:(54)Wnumerical=limε⟶1W=Wo+W1+W2.The convergence condition [43] due to the perturbation method is set as follows:(55)Wi+1Wi<1,wherei=0,1.


### 3.2.2. Iterative Quadrature

We solved the following iterative system:(56)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55Wi+1,xz+e1+e5Φ,xz=−ρω2U,i=0,1,2,(57)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44Wi+1,yz+e2+e4Φ,yz=−ρω2V,i=0,1,2,(58)C55+C13U,xz+C13+C44V,yz+C55−k2Wi+1,xx+C44−k2Wi+1,yy+C33Wi+1,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wi+1+k3Wi2Wi+1=−ρω2Wi+1,i=0,1,2,(59)e1+e5U,xz+e2+e4V,yz+e5Wi+1,xx+e4Wi+1,yy+e3Wi+1,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0,i=0,1,2.To getWo, we solved the following eigenvalue problem:(60)C11∑l=−NxNxBilxUljk+C66∑m=−NyNyBjmyUimk+C55∑n=−NzNzBknzUijn+C12+C66∑l=−NxNxAilx∑m=−NyNyAjmyVlmk+C13+C55∑l=−NxNxAilx∑n=−NzNzAknzW0,ljn+e1+e5∑l=−NxNxAilx∑n=−NzNzAknzΦljn=−ρω2U,(61)C12+C66∑l=−NxNxAilx∑m=−NyNyAjmyUlmk+C66∑l=−NxNxBilxVljk+C11∑m=−NyNyBjmyVimk+C44∑n=−NzNzBknzVijn+C13+C44∑m=−NyNyAjmy∑n=−NzNzAknzW0,imn+e2+e4∑m=−NyNyAjmy∑n=−NzNzAknzΦimn=−ρω2V,(62)C55+C13∑l=−NxNxAilx∑n=−NzNzAknzUljn+C13+C44∑m=−NyNyAjmy∑n=−NzNzAknzVimn+C55−k2∑l=−NxNxBilxW0,ljk+C44−k2∑m=−NyNyBjmyW0,imk+C33∑n=−NzNzBknzWijn+e5∑l=−NxNxBilxΦljk+e4∑m=−NyNyBjmyΦimk+e3∑n=−NzNzBknzΦijn+k1W0,ijk=−ρω2W0,(63)e1+e5∑l=−NxNxAilx∑n=−NzNzAknzUljn+e2+e4∑m=−NyNyAjmy∑n=−NzNzAknzVimn+e5∑l=−NxNxBilxW0,ljk+e4∑m=−NyNyBjmyW0,imk+e3∑n=−NzNzBknzW0,ijn−η1∑l=−NxNxBilxΦljk−η2∑m=−NyNyBjmyΦimk−η3∑n=−NzNzBknzΦljn=0.The boundary conditions (21)–(27) can also be approximated using two DQMs as follows:(1)
Simply supported (S):(64)Wijk=Vijk=C11∑l=−NxNxAilxUljk+C12∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e1∑n=−NzNzAknzΦijn=0,atx=0,a,Wijk=Uijk=C12∑l=−NxNxAilxUljk+C11∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e2∑n=−NzNzAknzΦijn=0,aty=0,b.(2)
Clamped (C):(65)Wijk=Vijk=Uijk=0,atx=0,ay=0,b.(3)
Free surface (F):(66)C11∑l=−NxNxAilxUljk+C12∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e1∑n=−NzNzAknzΦijn=0,C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=−NxNxAilxΦljk=C66∑l=−NxNxAilxVljk+∑m=−NyNyAjmyUimk=0,atx=0,a,C12∑l=−NxNxAilxUljk+C11∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e2∑n=−NzNzAknzΦijn=0,C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+e4∑m=−NyNyAjmyΦimk=C66∑l=−NxNxAilxVljk+∑m=−NyNyAjmyUimk=0,aty=0,b,Mechanical and electrical boundary conditions at the lower and upper surfaces of the composite are as follows:(67)C13∑l=−NxNxAilxUljk+C13∑m=−NyNyAjmyVimk+C33∑n=−NzNzAknzWijn+e3∑n=−NzNzAknzΦijn=C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=−NxNxAilxΦljk=0C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+e4∑m=−NyNyAjmyΦimk=e1∑l=−NxNxAilxUljk+e2∑m=−NyNyAjmyVimk+e3∑n=−NzNzAknzWijn−η3∑n=−NzNzAknzΦijn=0,atz=0,C13∑l=−NxNxAilxUljk+C13∑m=−NyNyAjmyVimk+C33∑n=−NzNzAknzWijn+e3∑n=−NzNzAknzΦijn=C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=1NxAilxΦljk=0C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+∑m=−NyNyAjmyΦimk=Φijk=0,atz=h.The continuity conditions between the interfaces of layers can be assumed as follows:(68)Wijk1=Wijk2,Vijk1=Vijk2,Uijk1=Uijk2,Φijk1=Φijk2.We have solved the generalized eigenvalue problem [65]:(69)KX=ω2MX,where K is the coefficient matrix of the previous system, M is the mass matrix and can be diagonal with 0 diagonal elements, and ω is the free vibration frequencies squared.For nontrivial solutions for (70), the following determinant should be 0:(70)K−ω2M=0.Equation (71) gives the natural frequencies for the composite plate.


## 3.2.1. Perturbation Method

Assuming the solution of equations (17)–(27) to be a power series that can be solved by a perturbation method of the second order [43] by substituting W=Wo+εW1+ε2W2+⋯ into (17)–(27) and equating the terms with the identical powers of ε, where ε is a perturbation parameter:(38)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55Wo,xz+εW1,xz+ε2W2,xz+e1+e5Φ,xz=−ρω2U,(39)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44Wo,yz+εW1,yz+ε2W2,yz+e2+e4Φ,yz=−ρω2V,(40)C55+C13U,xz+C13+C44V,yz+C55−k2Wo,xx+εW1,xx+ε2W2,xx+C44−k2Wo,yy+εW1,yy+ε2W2,yy+C33Wo,zz+εW1,zz+ε2W2,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wo+εW1+ε2W2+k3εWo+εW1+ε2W23=−ρω2Wo+εW1+ε2W2,(41)e1+e5U,xz+e2+e4V,yz+e5Wo,xx+εW1,xx+ε2W2,xx+e4Wo,yy+εW1,yy+ε2W2,yy+e3Wo,zz+εW1,zz+ε2W2,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.Equate the terms ofε0in (39)–(42):(42)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W0,xz+e1+e5Φ,xz=−ρω2U,(43)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W0,yz+e2+e4Φ,yz=−ρω2V,(44)C55+C13U,xz+C13+C44V,yz+C55−k2Wo,xx+C44−k2Wo,yy+C33Wo,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wo=−ρω2Wo,(45)e1+e5U,xz+e2+e4V,yz+e5W0,xx+e4W0,yy+e3W0,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.The previous system is an eigenvalue problem, which is solved to obtain the natural frequencies and the bending deflectionW0.Equate the terms ofε1 in (39)–(42):(46)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W1,xz+e1+e5Φ,xz=−ρω2U,(47)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W1,yz+e2+e4Φ,yz=−ρω2V,(48)C55+C13U,xz+C13+C44V,yz+C55−k2W1,xx+C44−k2W1,yy+C33W1,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1+ρω2W1=−k3Wo3,(49)e1+e5U,xz+e2+e4V,yz+e5W1,xx+e4W1,yy+e3W1,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0,The previous system is a bending problem, which is solved to obtainW1.Equate the terms ofε2 in (39)–(42):(50)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55W2,xz+e1+e5Φ,xz=−ρω2U,(51)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44W2,yz+e2+e4Φ,yz=−ρω2V,(52)C55+C13U,xz+C13+C44V,yz+C55−k2W2,xx+C44−k2W2,yy+C33W2,zz+e5Φ,xx+e4Φ,yy+e3Φ,zzk1+ρω2W2=−3k3Wo2W1,(53)e1+e5U,xz+e2+e4V,yz+e5W2,xx+e4W2,yy+e3W2,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0.The previous system is a bending problem, which is solved to obtainW2.Finally, the series solution can be written as follows:(54)Wnumerical=limε⟶1W=Wo+W1+W2.The convergence condition [43] due to the perturbation method is set as follows:(55)Wi+1Wi<1,wherei=0,1.


## 3.2.2. Iterative Quadrature

We solved the following iterative system:(56)C11U,xx+C66U,yy+C55U,zz+C12+C66V,xy+C13+C55Wi+1,xz+e1+e5Φ,xz=−ρω2U,i=0,1,2,(57)C66+C12U,xy+C66V,xx+C11V,yy+C44V,zz+C13+C44Wi+1,yz+e2+e4Φ,yz=−ρω2V,i=0,1,2,(58)C55+C13U,xz+C13+C44V,yz+C55−k2Wi+1,xx+C44−k2Wi+1,yy+C33Wi+1,zz+e5Φ,xx+e4Φ,yy+e3Φ,zz+k1Wi+1+k3Wi2Wi+1=−ρω2Wi+1,i=0,1,2,(59)e1+e5U,xz+e2+e4V,yz+e5Wi+1,xx+e4Wi+1,yy+e3Wi+1,zz−η1Φ,xx−η2Φ,yy−η3Φ,zz=0,i=0,1,2.To getWo, we solved the following eigenvalue problem:(60)C11∑l=−NxNxBilxUljk+C66∑m=−NyNyBjmyUimk+C55∑n=−NzNzBknzUijn+C12+C66∑l=−NxNxAilx∑m=−NyNyAjmyVlmk+C13+C55∑l=−NxNxAilx∑n=−NzNzAknzW0,ljn+e1+e5∑l=−NxNxAilx∑n=−NzNzAknzΦljn=−ρω2U,(61)C12+C66∑l=−NxNxAilx∑m=−NyNyAjmyUlmk+C66∑l=−NxNxBilxVljk+C11∑m=−NyNyBjmyVimk+C44∑n=−NzNzBknzVijn+C13+C44∑m=−NyNyAjmy∑n=−NzNzAknzW0,imn+e2+e4∑m=−NyNyAjmy∑n=−NzNzAknzΦimn=−ρω2V,(62)C55+C13∑l=−NxNxAilx∑n=−NzNzAknzUljn+C13+C44∑m=−NyNyAjmy∑n=−NzNzAknzVimn+C55−k2∑l=−NxNxBilxW0,ljk+C44−k2∑m=−NyNyBjmyW0,imk+C33∑n=−NzNzBknzWijn+e5∑l=−NxNxBilxΦljk+e4∑m=−NyNyBjmyΦimk+e3∑n=−NzNzBknzΦijn+k1W0,ijk=−ρω2W0,(63)e1+e5∑l=−NxNxAilx∑n=−NzNzAknzUljn+e2+e4∑m=−NyNyAjmy∑n=−NzNzAknzVimn+e5∑l=−NxNxBilxW0,ljk+e4∑m=−NyNyBjmyW0,imk+e3∑n=−NzNzBknzW0,ijn−η1∑l=−NxNxBilxΦljk−η2∑m=−NyNyBjmyΦimk−η3∑n=−NzNzBknzΦljn=0.The boundary conditions (21)–(27) can also be approximated using two DQMs as follows:(1)
Simply supported (S):(64)Wijk=Vijk=C11∑l=−NxNxAilxUljk+C12∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e1∑n=−NzNzAknzΦijn=0,atx=0,a,Wijk=Uijk=C12∑l=−NxNxAilxUljk+C11∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e2∑n=−NzNzAknzΦijn=0,aty=0,b.(2)
Clamped (C):(65)Wijk=Vijk=Uijk=0,atx=0,ay=0,b.(3)
Free surface (F):(66)C11∑l=−NxNxAilxUljk+C12∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e1∑n=−NzNzAknzΦijn=0,C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=−NxNxAilxΦljk=C66∑l=−NxNxAilxVljk+∑m=−NyNyAjmyUimk=0,atx=0,a,C12∑l=−NxNxAilxUljk+C11∑m=−NyNyAjmyVimk+C13∑n=−NzNzAknzWijn+e2∑n=−NzNzAknzΦijn=0,C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+e4∑m=−NyNyAjmyΦimk=C66∑l=−NxNxAilxVljk+∑m=−NyNyAjmyUimk=0,aty=0,b,Mechanical and electrical boundary conditions at the lower and upper surfaces of the composite are as follows:(67)C13∑l=−NxNxAilxUljk+C13∑m=−NyNyAjmyVimk+C33∑n=−NzNzAknzWijn+e3∑n=−NzNzAknzΦijn=C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=−NxNxAilxΦljk=0C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+e4∑m=−NyNyAjmyΦimk=e1∑l=−NxNxAilxUljk+e2∑m=−NyNyAjmyVimk+e3∑n=−NzNzAknzWijn−η3∑n=−NzNzAknzΦijn=0,atz=0,C13∑l=−NxNxAilxUljk+C13∑m=−NyNyAjmyVimk+C33∑n=−NzNzAknzWijn+e3∑n=−NzNzAknzΦijn=C55∑l=−NxNxAilxWljk+∑n=−NzNzAknzUijn+e5∑l=1NxAilxΦljk=0C44∑m=−NyNyAjmyWimk+∑n=−NzNzAknzVijn+∑m=−NyNyAjmyΦimk=Φijk=0,atz=h.The continuity conditions between the interfaces of layers can be assumed as follows:(68)Wijk1=Wijk2,Vijk1=Vijk2,Uijk1=Uijk2,Φijk1=Φijk2.We have solved the generalized eigenvalue problem [65]:(69)KX=ω2MX,where K is the coefficient matrix of the previous system, M is the mass matrix and can be diagonal with 0 diagonal elements, and ω is the free vibration frequencies squared.For nontrivial solutions for (70), the following determinant should be 0:(70)K−ω2M=0.Equation (71) gives the natural frequencies for the composite plate.


## 4. Numerical Results

The present numerical results demonstrate the convergence and efficiency of each one of the proposed schemes for vibration analysis of the piezoelectric composite plate Al/Al2O3 resting on a nonlinear elastic foundation. For all results, the boundary conditions (65)–(69) are augmented in the governing equations (39)–(42).For practical purposes, the field quantities are normalized as follows:(71)U∗,V∗,W∗=U,V,Wh,σx∗,σy∗,σz∗,τxy∗,τxz∗,τyz∗=σx,σy,σz,τxy,τxz,τyzE2,(72)Cij∗=CijE2,x∗=xa,y∗=yb,z∗=zh,h=hs+hc+ha,ha=hs=hp.(73)D=Ech3121−ν2,K1¯=k1a4D,K2¯=k2a2D,K3¯=k3a4h2D,where U∗,V∗,andW∗ are the normalized amplitudes of displacements, σx∗,σy∗,andσz∗ are the normalized amplitudes of stresses, τxy∗,τxz∗,andτyz∗ are the normalized amplitudes of shear stresses, haandhs are the thicknesses of the actuator and sensor, and D,K1,¯K2¯,andK3¯ are the flexural rigidity for the bottom layer of composite and the equivalent dimensionless elastic foundation constants.The computational characteristics of each scheme are adapted to reach accurate results with error of order ≤10−8. The obtained dimensionless frequencies Ω, ω¯, and β are evaluated as follows: Ω=ωhρ/E,ω¯=ωb/π2ρh/D,andβ=ωa2ρh/D, where ρ and E are the density and Young’s modulus, respectively.For the present results, material parameters for the composite are listed in Table1.Table 1
Material property of composite piezoelectric plate [59, 61].
Material propertyYoung’s moduli (GPa)Shear moduli (GPa)Poisson’s ratiosDensity (kg/m3)E2E1E3G12G13G23ν12ν13ν23ρ725E2E20.5E2G120.2E20.250.030.41600Metal (aluminum, Al)70———0.32702Ceramic (alumina, Al2O3)380———0.33800Effective elastic (GPa)C11C12C13C22C23C33C44C55C66Sensor PZT-413978741397411525.625.630.5Sensor BaTiO3166777816678162434344.5Actuator Ba2NaNb5O15239104524752135656676Actuator PZT-5A121777712111121212123Material propertyPiezoelectric constants (C/m2)Dielectric constants (F/m) ∗ 10−9Density (kg/m3)e1e2e3e4e5η1η2η3ρSensor PZT-4−5.2−5.215.112.712.76.56.55.67500Sensor BaTiO3−4.4−4.418.611.611.611.211.212.65700Actuator Ba2NaNb5O15−0.4−0.34.33.42.819.62.010.285300Actuator PZT-5A−5.4−5.415.83.43.48.118.117.342330For the sinc DQ scheme, the problem is solved over a regular grid, ranging from 3 × 5 × 5 to 11 × 5 × 5. Table2 shows the convergence of the obtained results for the isotropic plate, which is in agreement with the exact results in [68] over grid size ≥7 × 5 × 5.Table 2
Comparison between the normalized fundamental frequenciesβ=ωa2ρmh/Dm using the sinc DQM, grid size N, and the previous exact and numerical frequencies for the simply supported square isotropic plate (a/h = 5, hc/hp=200, K1¯=1000, K2¯=100, and K3¯=0).
ResultsNormalized frequenciesβ1β2β3β4β5β6NSinc DQM355.9218103.4553124.1464161.3903177.5293186.4058555.9721103.5484124.2581161.5355177.6891186.5736756.0369103.6682124.2982161.5391177.6930186.5777956.0369103.6683124.2982161.5391177.6930186.57771156.0369103.6683124.2982161.5391177.6930186.5777Analytical [68]56.0359—————RHSDT [69]56.0310—————NSDT [60]56.0311—————Execution time (sec)5.615225 overN = 7 × 5 × 5RHSDT: refined hyperbolic shear deformation theory; NSDT: new shear deformation theory; DQM: differential quadrature method.For DSCDQ scheme based on Delta Lagrange kernel, the problem is also solved over a uniform grid ranging from 3 × 5 × 5 to 11 × 5 × 5. The bandwidth 2M + 1 ranges from 3 to 11.Table3 shows the convergence of the obtained fundamental frequency, which is in agreement with the exact results in [68] over grid size ≥7 × 5 × 5 and bandwidth ≥5. Table 4 shows that the obtained results are more accurate than those obtained using approximated theories [60, 69]. This table also shows that the execution time of DSCDQM-DLK was less than that of sinc DQM.Table 3
Comparison between the normalized fundamental frequency using DSCDQM-DLK, bandwidth (2M + 1), and grid size N for the simply supported square isotropic plate (hc/hp=200, a/h=5, K1¯=1000, K2¯=100, and K3¯=0).
Fundamental frequencyDSCDQM-DLKBandwidthN3579112M + 1 = 355.984456.084956.068556.068556.06852M + 1 = 555.952456.034956.032056.032056.03202M + 1 = 756.179856.065056.032056.032056.03202M + 1 = 955.938456.287756.032056.032056.03202M + 1 = 1155.938456.287756.032056.032056.0320Table 4
Comparison between the normalized frequenciesβ=ωa2ρmh/Dm, using DSCDQM-DLK, grid sizes N, and the previous exact and numerical frequencies for the simply supported square isotropic plate (2M + 1 = 5, hc/hp=200, a/h=5, K1¯=1000, K2¯=100, and K3¯=0).
Normalized frequenciesβ1β2β3β4β5β6ResultsNDSCDQM-DLK355.9524103.5119124.2142161.4660177.6126186.4932556.0349103.6085124.3302161.5049177.6553186.5380756.0320103.6031124.3237161.4964177.6444186.5266956.0320103.6031124.3237161.4964177.6444186.5266Analytical [68]56.0359—————RHSDT [69]56.0310—————NSDT [60]56.0311—————Execution time (sec)5.611001 overN = 7 × 5 × 5For the DSCDQ scheme based on the regularized Shannon kernel (RSK), the problem was also solved over a uniform grid ranging from 3 × 5 × 5 to 9 × 5 × 5. The bandwidth 2M + 1 ranged from 3 to 11 and the regularization parameter σ = r hx ranged from 1.0 hx to 2.3 hx, where hx = 1/N − 1. Table 5 shows the convergence of the obtained fundamental frequency to the exact results in [68] over grid size ≥5 × 5 × 5, bandwidth ≥3, and regularization parameter σ = 1.9 hx.Table 5
Comparison between the normalized fundamental frequencyβ using the discrete singular convolution differential quadrature method on regularized Shannon kernel (DSCDQM-RSK), bandwidth (2M + 1), regularization parameter σ, and grid size N for the simply supported square isotropic plate (hc/hp=200, a/h=5, K1¯=1000, K2¯=100, and K3¯=0).
Fundamental frequencyRegularization parameterDSCDQM-RSKN2M + 1σ = 1 × hxσ = 1.5 × hxσ = 1.9 × hxσ = 2.3 × hx3355.956455.946855.985356.0821555.956155.899155.978056.0706755.956155.899155.978056.07065355.875855.971556.031856.0318555.875855.971556.031856.0318755.875855.971556.031856.03187355.850755.939256.031856.0318555.956355.979056.031856.0318755.992655.991756.031856.03189355.925955.985756.031856.0318555.982355.973356.031856.0318755.950656.019356.031856.0318Table6 shows that the obtained results from DQ schemes were more accurate than other models in the literature [60, 69]. Further, the execution time of this scheme is the least. Therefore, the DSCDQM-RSK scheme is the best choice among the examined quadrature schemes for vibration analysis of piezoelectric composite materials resting on nonlinear elastic foundations. Furthermore, a parametric study is introduced to investigate the influence of foundation parameters, elastic and geometric characteristics of the composite, and types of material on the values of natural frequencies and mode shapes. Tables 7 and 8 show the natural frequency for isotropic plates. These tables show the effect of elastic foundation parameters on the frequency for different boundary conditions. The natural frequency increased with the increase in foundation parameters K1¯ and K2¯. Additionally, the SCSC plate had the highest frequency and the SSSF had the lowest frequency. Moreover, the efficiency of the nonlinear treatment methods was investigated. The execution time for the iterative quadrature method was less than the perturbation method, with a factor close to 1.5 until 1.9 for all boundary conditions used in Table 8.Table 6
Comparison between the normalized frequenciesβ=ωa2ρmh/Dm using DSCDQM-RSK, grid size, and the previous exact and numerical frequencies for the simply supported square isotropic plate (2M + 1 = 3, σ = 1.9 ∗ hx, hc/hp=200, a/h=5, K1¯=1000, K2¯=100, and K3¯=0).
Normalized frequenciesβ1β2β3β4β5β6ResultsNDSCDQM-RSK355.9853103.5168124.2201161.4861177.6347186.5164556.0318103.6027124.3232161.6201177.7821186.6712756.0318103.6027124.3232161.6201177.7821186.6712956.0318103.6027124.3232161.6201177.7821186.6712Analytical [68]56.0359—————RHSDT [69]56.0310—————NSDT [60]56.0311—————Execution time (sec)1.933552 overN = 5 × 5 × 5Table 7
Comparison between the normalized fundamental frequencyβ=ωa2ρmh/Dm and the elastic foundation parameter for the simply supported squared plate (hc/hp=200 and a/h = 5).
K1¯, K2¯Nonlinear treatment methodsK3¯ = 0K3¯ = 50K3¯ = 1000, 0Perturbation17.459017.459017.4590Execution time (sec)0.4807631.23461.2210Iterative quadrature17.459017.459017.4590Execution time (sec)0.4807630.7434080.818601Analytical [68]17.5055——100, 10Perturbation24.280724.280724.2807Execution time (sec)0.4699891.17441.2081Iterative quadrature24.280724.280724.2807Execution time (sec)0.4699890.6291480.651756Analytical [68]24.3074——1000, 100Perturbation56.031856.031856.0318Execution time (sec)0.4394781.14961.1351Iterative quadrature56.031856.031856.0318Execution time (sec)0.4394780.6204710.572076Analytical [68]56.0359——Table 8
Comparison between the normalized fundamental frequencyβ, elastic foundation parameter, and different boundary conditions for the square plate (hc/hp=200, a/h=5, and K3¯ = 50).
K1¯K2¯Nonlinear treatment methodsSSSSSCSCSSSCSSSF00Perturbation19.734928.915723.639811.6883Iterative quadrature19.734928.915723.639811.6883100Perturbation48.632954.650551.336537.1656Iterative quadrature48.632954.650551.336537.16561000Perturbation141.9431146.7559144.2712111.7952Execution time (sec)1.241521.21581.221581.187421Iterative quadrature141.9425146.7553144.2706111.7947Execution time (sec)0.7845090.703550.764870.694571000Perturbation22.140230.674125.681815.3870Iterative quadrature22.140230.674125.681815.3870100Perturbation49.632455.568952.292538.4718Iterative quadrature49.632455.568952.292538.47181000Perturbation142.2904147.1665144.7015112.3010Execution time (sec)1.241151.11871.240051.12571Iterative quadrature142.3024147.1790144.7098112.3075Execution time (sec)0.777850.702540.757140.7047810000Perturbation37.279142.829339.488933.7090Iterative quadrature37.279142.829339.488933.7090100Perturbation57.956963.148560.854748.6837Iterative quadrature57.956963.148560.854748.68371000Perturbation145.3461150.1200147.6716116.1315Execution time (sec)1.142871.112761.014510.99475Iterative quadrature145.4915150.2702147.6863116.2471Execution time (sec)0.754150.701250.745830.69573Tables8–10 and Figures 2–5 show that the natural frequencies increase with increasing side-to-thickness ratio (a/h), foundation parameters K1¯,K2¯,K3¯, total thickness, Young’s modulus gradation ratio (E1/E2), shear modulus gradation ratio (G13/G12), and number of layers. The natural frequencies converge at a total thickness (h = 0.075) in the CCCC plate and a total thickness (h = 0.035) in the SSSSS plate, as shown in Tables 10 and 11. Figure 6 shows the natural frequencies decrease with the decrease of the piezoelectric layer thickness (hc/hp).Table 9
Comparison between the normalized frequenciesβ=ωa2ρch/Dc, thickness ratio (a/h), elastic foundation parameter, and different boundary conditions for the square plate (hc/hp=200, K2¯=10, and K3¯ = 50).
a/hK1¯SCSCSCSSω1¯ω2¯ω3¯ω1¯ω2¯ω3¯10/302.31503.42763.67422.17622.45723.5834102.33283.42803.63242.19052.45743.59091002.48343.43113.63252.34152.46163.712410003.42983.43923.63263.32213.12123.8330201.76292.23832.72141.70721.81672.6100101.77522.28832.76231.71321.82322.61121001.85892.28852.76281.79211.86722.619410002.09772.33842.74562.04921.96882.6414502.38974.71397.75062.24563.39407.5180102.41164.72407.75622.33923.92097.52471002.60014.81357.80542.41803.99897.591810003.97085.59448.25942.41814.00157.5921Table 10
Comparison between the normalized frequencies and the thicknessh for the clamped plate (a/h=4,k1¯=1000,k2¯=100,andk3¯=50).
Normalized frequenciesβ1β2β3β4β5β6H0.0579.648164.334211.562237.039301.984327.4230.0679.918164.853211.828237.510302.594327.9960.0779.918164.854211.853237.512302.594328.0100.07579.918164.854211.853237.512302.594328.0100.0879.918164.854211.853237.512302.594328.0100.0979.918164.854211.853237.512302.594328.0100.179.918164.854211.853237.512302.594328.010Figure 2
Variation of the fundamental frequency (Ω) with Young’s modulus (E1/E2), different boundary conditions, and different materials for the square plate (K1¯=1000, K2¯=100, K3¯=50, hc/hp=45, and a/h=10). (a) Sensor is PZT-4 and actuator is Ba2NaNb5O15; (b) sensor is BaTiO3 and actuator is PZT-5A.
(a)(b)Figure 3
Variation of the fundamental frequency (Ω) with shear modulus (G13/G12), different boundary conditions, and different materials of the square plate (K1¯=1000,K2¯=100,K3¯=50, hc/hp=45,anda/h=10). (a) Sensor is PZT-4 and actuator is Ba2NaNb5O15; (b) sensor is BaTiO3 and actuator is PZT-5A.
(a)(b)Figure 4
Variation of the fundamental frequency (Ω) with the number of layers, different boundary conditions, and different materials of the square plate (K1¯=1000,K2¯=100,K3¯=100, hc/hp=45,anda/h=5). (a) Sensor is PZT-4 and actuator is Ba2NaNb5O15; (b) sensor is BaTiO3 and actuator is PZT-5A.
(a)(b)Figure 5
Variation of the fundamental frequency (Ω) with Young's modulus (E1/E2), shear modulus (G13/G12), and different number of layers for the different materials (sensor is BaTiO3 and actuator is PZT-5A) of the square plate (K1¯=1000,K2¯=100,K3¯=100, hc/hp=45,anda/h=10).
(a)(b)Table 11
Comparison between the normalized frequencies and thicknessh for the simply supported plate (a/h=4,k1¯=1000,k2¯=100,andk3¯=50).
Normalized frequenciesβ1β2β3β4β5β6H0.0151.612292.90196113.5468144.5142165.159180.64270.0251.787493.21732113.9323145.0047165.7197181.25590.0351.787593.2175113.9325145.005165.7200181.25630.03551.787593.2175113.9325145.005165.7200181.25630.0451.787593.2175113.9325145.005165.7200181.25630.0551.787593.2175113.9325145.005165.7200181.25630.0751.787593.2175113.9325145.005165.7200181.2563Figure 6
Variation of the fundamental frequency (Ω) with the composite layer thickness to piezoelectric thickness ratio (hc/hp), length-to-thickness ratio (a/h), and different boundary conditions for the square plate (K1¯=1000,K2¯=100,K3¯=50, and E1/E2=25). (a) a/h = 5 and (b) a/h = 10.
(a)(b)Furthermore, Figure7 shows that the natural frequencies decreased with an increase in the aspect ratio (a/b) at different values (a/h). Moreover, the natural frequencies remained constant when hc/hp > 20 in all different edge conditions. Figures 2, 3, and 8 show the effect of the material type on the natural frequency. It can be seen that the sensor was more significant than the actuators. The natural frequency was almost unchanged when a/h > 40. Additionally, the influence of the BaTiO3 material was more affected by the natural frequencies than the PZT-4 material. Furthermore, Figures 9–12 show the first three mode shapes of normalized transverse displacement (w) and electric potential (φ). From the previous figures, it can be seen that the normalized transverse displacement was at a maximum for the CCCF plate and was at a minimum for the CCCC plate. Additionally, the electric potential was at a maximum for the CCCF plate and at a minimum for the CSCS plate. Furthermore, the previous figures show that the amplitude of displacement and electric potential increased with the increase of the nonlinear Winkler parameter K3¯.Figure 7
Variation of the fundamental frequency (Ω) with the aspect ratio (a/b), different values of a/h, and different boundary conditions for the square plate (K1¯=1000,K2¯=100,K3¯=100, hc/hp=25,andE1/E2=25). (a) SSSS, (b) CSCS, (c) SCSC, and (d) CCCF.
(a)(b)(c)(d)Figure 8
Variation of the fundamental frequency (Ω) with the side-to-thickness ratio (a/h), different boundary conditions, and different materials for the square plate (K1¯=1000, K2¯=100,K3¯=100, and E1/E2=25). (a) Sensor is BaTiO3 and actuator is Ba2NaNb5O15; (b) sensor is BaTiO3 and actuator is PZT-5A.
(a)(b)Figure 9
Variation of normalized mode shape deflectionW with the time for the first three modes at different boundary conditions for the square plate (K1¯=1000,K2¯=100,K3¯=0, hc/hp=25,E1/E2=25,anda/h=10). (a) SSSS, (b) CCCC, (c) SCSC, and (d) CCCF.
(a)(b)(c)(d)Figure 10
Variation of normalized mode shape deflectionW with the time for the first three modes at different boundary conditions for the square plate (K1¯=1000,K2¯=100,K3¯=100, hc/hp=25,E1/E2=25,anda/h=10). (a) SSSS, (b) CCCC, (c) SCSC, and (d) CCCF.
(a)(b)(c)(d)Figure 11
Variation of electric potential with the time for the first three modes at different boundary conditions for the square plate (K1¯=1000, K2¯=100, K3¯=0, hc/hp=25, E1/E2=25, and a/h=10). (a) SSSS, (b) CCCC, (c) SCSC, and (d) CCCF.
(a)(b)(c)(d)Figure 12
Variation of electric potential with the time for the first three modes at different boundary conditions for the square plate (K1¯=1000,K2¯=100,K3¯=100, hc/hp=25,E1/E2=25,anda/h=10). (a) SSSS, (b) CCCC, (c) SCSC, and (d) CCCF.
(a)(b)(c)(d)


## 5. Conclusions

Accurate quadrature schemes have been applied for free vibration analysis of piezoelectric composite plates with and without nonlinear elastic foundations for different types of boundary conditions. Numerical results were given to demonstrate the convergence, accuracy, and efficiency of the present methods. A MATLAB program was designed for each one, such that the maximum error (comparing with the previous exact results) was ≤10−8. Additionally, the execution time for each scheme was determined. It is concluded that the convergence speed for the discrete singular convolution differential quadrature method based on regularized Shannon kernel (DSCDQM-RSK) with grid size ≥5 × 5 × 5, bandwidth 2M + 1 ≥ 3, and regularization parameter σ = 1.9 hx is very fast and has the best accurate efficient results for the concerned problem (three-layered piezoelectric composite with total thickness = 0.075 in the CCCC plate and h = 0.035 in the SSSS plate). The efficiency of the iterative quadrature method is better than the perturbation technique. The influences of linear and nonlinear elastic foundation on the natural frequency and higher modes of frequency are investigated. The effect of the Pasternak parameter is larger than the effect of the linear and nonlinear Winkler coefficients on the natural frequencies. The natural frequency does not change when varying the nonlinear Winkler parameter. Additionally, the natural frequency converges at a larger value of the linear Winkler parameter. A parametric study is also introduced to investigate the influence of elastic and geometric characteristics of the composite and types of materials of the vibrated plate on the results. The nondimensional frequency for the material type BaTiO3 is larger than that for the PZT-4 material. Further, the natural frequency decreases with the decrease in the thickness of the piezo material or with the increase in the composite thickness until the value of ratio hc/hp closes to 10, so the natural frequency remains constant. The thinner composite has frequencies larger than thick composite. The composite plate with CCCC conditions has the minimum normalized transverse displacement, whereas it has the maximum value for CCCF conditions. Additionally, the electric potential is the maximum for the CCCF plate and has a minimum value for the CSCS plate. Moreover, the vertical displacement increases slightly when enlarging the nonlinear parameter. Finally, the aim is for these results to be useful for modern engineering smart structures.



---

*Source: 7592302-2020-07-15.xml*