Three-dimensional static analysis of a viscoelastic rectangular functionally graded material plate embedded between piezoelectric sensor and actuator layers

Abstract

The three-dimensional bending behavior of a viscoelastic functionally graded material (FGM) layer embedded between piezoelectric layers and subjected to an electric field as well as a uniform transverse pressure is studied. An analytical solution is computed for a simply supported viscoelastic smart FGM plate using the state space technique along the thickness direction and a Fourier expansion along the in-plane coordinates. The governing differential equations in the time domain are transformed to the Laplace domain, solved in the Laplace domain, and transformed back to the time domain using the inverse Laplace transform. In the present study, the relaxation modulus of the FGM layer is assumed in the form of a Prony series, and it varies according to the power law in the thickness direction. The validity of the proposed approach is assessed by comparison of the numerical results of the approach with those of published works in the literature. The effects of geometric dimensions, electromechanical loads, and relaxation time constant on the behavior of the smart plate are investigated.

Keywords:

• Functionally graded material
• static
• viscoelastic
• piezoelectric
• plates

Introduction

Functionally graded materials (FGMs) represent specific types of composite materials whose properties vary gradually along a specific spatial direction as a result of their non-uniform microstructures with continuously graded macro-properties (Apalak and Demirbas 2019). The resistance of FGMs in high temperatures as a thermal barrier makes them a desirable option in environments such as rocket engines and space bodies (Marin 2005).

Brischetto et al. (2008) used Extending Unified Formulation introduced by Carrera to investigate deformations of a simply supported, functionally graded, rectangular plate subjected to thermo- mechanical loadings using Fourier series solution. Based on Carrera Unified Formulation (CUF), Bharati et al. (2020) analyzed flutter of a laminated composite structures under aerodynamics loading conditions using the principle of virtual displacement along with finite element method. Stress distributions in the geometrical nonlinear regime of functionally graded structures was evaluated by Demirbas et al. (2020) in the framework of CUF employing Newton-Raphson technique as well as Finite Element Method. Based on higher order shear deformation theory, Sahoo et al. (2020) developed an eigenvalue-based approach to predict the critical buckling temperature of graded sandwich structures. Sahoo et al. (2021) studied numerically thermal eigenvalue responses of a functionally graded sandwich shell structure. Effect of grading pattern and porosity on the eigen characteristics of porous functionally graded structure was studied by Ramteke et al. (2019) using isoparametric finite element technique. Ramteke et al. (2020) carried out static analysis of a two directional graded structure using a commercial FE method. Bending behavior of a functionally graded structure with variable grading patterns as well as variable porosity was studied by Ramteke et al. (2021) using finite element method.

At elevated temperatures, some materials, including the FGMs, tend to have creep and relaxation behaviors, which necessitates the incorporation of the theory of viscoelasticity, in which the time-dependent response cannot be considered negligible because of their viscoelastic behavior (Lakes and Lakes 2009). Widespread applications of FGMs in high-temperature environments causes a viscoelastic constitutive model to be more accurate compared with an elastic model when it comes to the analysis of such materials used in practical applications. Exhibiting electromechanical coupling effects between the electric field and the elastic deformation, piezoelectric materials combined with FGMs represent a smart or intelligent system. Piezoelectric layers, used as actuator and sensor layers, enable the system to respond to changing operating environments. Liew et al. (2001) investigated the three-dimensional transient response in FG rectangular plates equipped with piezoelectric layers. Phung-Van et al. (2013) used a Mindlin plate theory to study the static and dynamic behavior of laminated plates with piezoelectric layers. Khalili, Malekzadeh, and Mittal (2005) analyzed the static and dynamic behavior of a simply supported plate using the Fourier series technique and Stokes’s transformation. Vel and Batra (2000) investigated the bending behavior of rectangular laminated piezoelectric plates analytically using the Eshelby-Stroh formalism. Mitchell and Reddy (1995) developed a refined theory of laminated composite plates with piezoelectric laminae within the framework of linear piezoelectricity, in which they included the coupling between mechanical deformations and the charge equations of electrostatics. Fernandes and Pouget (2001) presented an efficient two-dimensional approach for piezoelectric plates in the framework of the linear theory of piezoelectricity using a variational method. Based on higher-order shear deformation theory and using the principle of virtual displacements, Khetib et al. (2019) studied the vibration behavior of an FGM plate resting on a viscoelastic Winkler-Pasternak foundation. Based on concepts of Boltzmann and Volterra, Sofiyev, Zerin, and Kuruoglu (2020) investigated free vibrations and the dynamic stability of an FGM viscoelastic plate under a compressive load and resting on elastic foundations using the Galerkin and Laplace integral transform. Buckling and free vibration analysis of sandwich plates with a viscoelastic core and FGM face sheets under high-temperature environment were carried out by Joseph and Mohanty (2019) in the framework of first-order shear deformation theory using the finite element method. Shokrgozar et al. (2020) studied, analytically and numerically, bending behavior of cylindrical microshell rested on viscoelastic foundation using Fourier series state-space technique and generalized differential quadrature method (GDQM).

Despite the fact that the analytical solutions are well accepted in the engineering domain, for cases such as viscoelastic structures in which the governing equations are hard to be handled analytically due to the structure shape and complex boundary conditions, the numerical solutions would be adopted (Sahoo et al. 2020). Norouzi and Alibeigloo (2016, 2018) investigated the three-dimensional bending behavior of a viscoelastic FGM cylindrical panel subjected to thermal or mechanical loads using a combined state space and Laplace transform approach. Active vibration control of FGM plates with integrated piezoelectric sensors and actuators in a thermal environment was studied using the finite element method in He et al. (2001) and Liew et al. (2003). Jermsittiparsert et al. (2020) used GDQ to investigate buckling, critical voltage, and vibration response of a thermally affected graphene nanoplatelet reinforced composite (GPLRC) micro disk coupled with a piezoelectric actuator.

Based on Reddy’s higher-order shear deformation plate theory, Yang, Kitipornchai, and Liew (2003) investigated the vibration behavior of FGM plates with integrated piezoelectric layers in a thermal environment using a one-dimensional differential quadrature and a Galerkin technique. Tekin and Kadıoğlu (2017) developed a mixed-type finite element model as well as a corresponding functional in the Laplace–Carson domain to analysis of quasi-static and dynamic behavior of viscoelastic Kirchhoff plates. They used Dubner and Abate numerical inverse Laplace to convert the solution to time domain. Paulino and Jin (2001) used the viscoelastic correspondence principle and the finite element method to simulate viscoelastic FGMs subjected to mechanical and thermal loads. Thermoelastic buckling and vibration analysis of an FGM sandwich beam with a constrained viscoelastic core was performed by Bhangale and Ganesan (2006) using the finite element method. To study the generalized theories of magneto-thermo-viscoelasticity, a rotating FGM plate was modeled numerically using the dual reciprocity boundary element method in Fahmy (2013).

Based on modified couple stress theory, vibration and stability analysis of a viscoelastic micro-shell conveying fluid in thermal environment was carried out by Rashvand, Alibeigloo, and Safarpour (2020). Safarpour, Rahimi, and Alibeigloo (2020) studied graphene reinforced composite truncated conical shells, cylindrical shells and annular plates sem-analytically using state space DQM. Three-dimensional bending and free vibration analysis of laminated composite plates regarding the effects of a viscoelastic interface was conducted by Alibeigloo (2015) using a state-space technique as well as Laplace transform.

The governing equations based on three-dimensional theory of elasticity are exact due to no simplifying assumptions in the extraction process of the equations. So, solutions of those equations can be used to assess the validity of results based on approximate conventional two-dimensional formulation. Moreover, in theory of elasticity, transverse normal and shear stresses which cause delamination and crack initiation at the interface of multi layered structures are determined exactly. Therefore, in this paper we use elasticity formulation to investigate, analytically, bending behavior of viscoelastic smart structures under electromechanical loads that has not yet been reported in literature.

The viscoelastic constitutive equations are expressed using the Boltzmann superposition principle, while the material properties of the FGM layer are assumed to vary in the thickness direction according to a power law. We use Fourier series approach along the in-plane coordinates and state-space technique across the thickness direction for space domain and Laplace transform for time domain.

Mathematical formulation

A simply supported sandwich rectangular plate with the geometry as shown in Fig. 1 is studied. The sandwich plate consists of three layers, namely, the top piezoelectric layer as the actuator, a viscoelastic FGM core layer in the middle, and the bottom piezoelectric layer as the sensor. $\mathit{a}$ and $\mathit{b}$ are the length and width of the plate along $\mathit{x}$ and $\mathit{y}$ directions, respectively, and the thickness of the total plate is $\mathit{h};$ the origin of the coordinate system is located on the bottom of the plate. The plate is subjected to the uniform pressure, ${\mathit{P}}_0,$ on the top surface as well as the electric voltage, $\mathit{V},$ at the top surface of the piezoelectric layer.

Figure 1. Sandwich plate with FGM layer sandwiched between two piezoelectric layers.

Three-dimensional governing equilibrium equations in the absence of body force are as follows: (1) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{fx}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{x}}{+\mathit{\sigma }}_{\mathit{\text{fy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{y}}{+\mathit{\sigma }}_{\mathit{\text{fz}},\mathit{z}}=0.\end{array}$$(1)

Most materials, which respond elastically at room temperature, tend to exhibit significant viscoelastic behavior at an elevated temperature. Since FGMs are usually used in high-temperature environments, one must consider viscoelasticity in order to capture the actual material behavior (Alibeigloo 2015). For linear viscoelastic materials, the three-dimensional constitutive relations can be represented by the Boltzmann superposition principle as Li et al. (2019) and Norouzi and Alibeigloo (2016), (2) $${\mathit{\sigma }}_{\mathit{i}}=\underset{-\infty }{\overset{\mathit{t}}{\int }}{\mathit{C}}_{\mathit{\text{ij}}}\left(\mathit{t}-\mathit{\tau }\right)\frac{\mathit{d}}{\mathit{\text{dτ}}}\left({\mathit{\epsilon }}_{\mathit{i}}\right)\mathit{d}\mathit{\tau },\hspace{1em}\mathit{i},\mathit{j}=1,2,\dots ,6,$$(2) where ${\mathit{\sigma }}_{\mathit{i}}(i=1,2,3)$ represents normal stresses, ${\mathit{\sigma }}_x,{\mathit{\sigma }}_y,$ and ${\mathit{\sigma }}_z;$${\mathit{\sigma }}_{\mathit{i}}(i=4,5,6)$ represents shear stresses, ${\tau }_{xy},{\tau }_{xz},$ and ${\tau }_{yz}.$ And similarly, ${\epsilon }_{\mathit{i}}(i=1,2,3)$ represents normal strains, ${\epsilon }_x,{\epsilon }_y,$ and ${\epsilon }_z;$ ${\epsilon }_{\mathit{i}}(i=4,5,6)$ represents shear strains, ${\gamma }_{xy},{\gamma }_{xz},$ and ${\gamma }_{yz}.$${\mathit{C}}_{\mathit{\text{ij}}}$ represents relaxation moduli coefficients, which depend on time in viscoelastic materials, defined in Eq. (A-1)(1) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{fx}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{x}}{+\mathit{\sigma }}_{\mathit{\text{fy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{y}}{+\mathit{\sigma }}_{\mathit{\text{fz}},\mathit{z}}=0.\end{array}$$(1) .

FGM layer

Formulation in Laplace domain

A set of integro-differential equations in the time domain is derived by substitution of Eq. (2)(2) $${\mathit{\sigma }}_{\mathit{i}}=\underset{-\infty }{\overset{\mathit{t}}{\int }}{\mathit{C}}_{\mathit{\text{ij}}}\left(\mathit{t}-\mathit{\tau }\right)\frac{\mathit{d}}{\mathit{\text{dτ}}}\left({\mathit{\epsilon }}_{\mathit{i}}\right)\mathit{d}\mathit{\tau },\hspace{1em}\mathit{i},\mathit{j}=1,2,\dots ,6,$$(2) into Eq. (1)(1) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{fx}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{x}}{+\mathit{\sigma }}_{\mathit{\text{fy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{y}}{+\mathit{\sigma }}_{\mathit{\text{fz}},\mathit{z}}=0.\end{array}$$(1) . To avoid the complexities accompanied with the solution of those equations, before the implementation, Eq. (2)(2) $${\mathit{\sigma }}_{\mathit{i}}=\underset{-\infty }{\overset{\mathit{t}}{\int }}{\mathit{C}}_{\mathit{\text{ij}}}\left(\mathit{t}-\mathit{\tau }\right)\frac{\mathit{d}}{\mathit{\text{dτ}}}\left({\mathit{\epsilon }}_{\mathit{i}}\right)\mathit{d}\mathit{\tau },\hspace{1em}\mathit{i},\mathit{j}=1,2,\dots ,6,$$(2) is transformed into the Laplace domain by using Laplace transform (Norouzi and Alibeigloo 2018): (3) $${\widehat{\sigma }}_{\mathit{i}}{\mathit{\text{(s)}}=\mathit{s}\widehat{C}}_{\mathit{\text{ij}}}\left(\mathit{s}\right){\widehat{\epsilon }}_{\mathit{i}}\mathit{\text{(s)}},$$(3) where the caret over each variable indicates the Laplace transformation of the variable.

Assuming a constant Poisson ratio, the linear viscoelastic constitutive relation in term of displacement components for FGM layer of the plate is expressed as (4) $$\begin{array}{c}{\widehat{\sigma }}_{fx}=\frac{s\widehat{E}\left(z,s\right)}{\left(1+\nu \right)\left(1-2\nu \right)}\left[\left(1-\nu \right){\widehat{u}}_{f,x}+\nu {\widehat{v}}_{f,y}+\nu {\widehat{w}}_{f,z}\right]\mathit{}\mathit{}\\ {\widehat{\sigma }}_{fy}=\frac{s\widehat{E}\left(z,s\right)}{\left(1+\nu \right)\left(1-2\nu \right)}\left[\nu {\widehat{u}}_{f,x}+\left(1-\nu \right){\widehat{v}}_{f,y}+\nu {\widehat{w}}_{f,z}\right]\\ {\widehat{\sigma }}_{fz}=\frac{s\widehat{E}\left(z,s\right)}{\left(1+\nu \right)\left(1-2\nu \right)}\left[\nu {\widehat{u}}_{f,x}+\nu {\widehat{v}}_{f,y}+\left(1-\nu \right){\widehat{w}}_{f,z}\right]\\ {\widehat{\tau }}_{\mathit{fxy}}=\frac{s\widehat{E}\left(z,s\right)}{2\left(1+\nu \right)}\left({\widehat{u}}_{f,y}+{\widehat{v}}_{f,x}\right)\\ {\widehat{\tau }}_{\mathit{fxz}}=\frac{s\widehat{E}\left(z,s\right)}{2\left(1+\nu \right)}\left({\widehat{u}}_{f,z}+{\widehat{w}}_{f,x}\right)\mathit{}\\ {\widehat{\tau }}_{\mathit{fyz}}=\frac{s\widehat{E}\left(z,s\right)}{2\left(1+\nu \right)}\left({\widehat{v}}_{f,z}+{\widehat{w}}_{f,y}\right),\end{array}$$(4) here $\mathit{f}$ indicates the FGM layer.

From Eqs. (1)(1) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{fx}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxy}},\mathit{x}}{+\mathit{\sigma }}_{\mathit{\text{fy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{fxz}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{fyz}},\mathit{y}}{+\mathit{\sigma }}_{\mathit{\text{fz}},\mathit{z}}=0.\end{array}$$(1) , (3)(3) $${\widehat{\sigma }}_{\mathit{i}}{\mathit{\text{(s)}}=\mathit{s}\widehat{C}}_{\mathit{\text{ij}}}\left(\mathit{s}\right){\widehat{\epsilon }}_{\mathit{i}}\mathit{\text{(s)}},$$(3) , and (4)(4) $$\begin{array}{c}{\widehat{\sigma }}_{fx}=\frac{s\widehat{E}\left(z,s\right)}{\left(1+\nu \right)\left(1-2\nu \right)}\left[\left(1-\nu \right){\widehat{u}}_{f,x}+\nu {\widehat{v}}_{f,y}+\nu {\widehat{w}}_{f,z}\right]\mathit{}\mathit{}\\ {\widehat{\sigma }}_{fy}=\frac{s\widehat{E}\left(z,s\right)}{\left(1+\nu \right)\left(1-2\nu \right)}\left[\nu {\widehat{u}}_{f,x}+\left(1-\nu \right){\widehat{v}}_{f,y}+\nu {\widehat{w}}_{f,z}\right]\\ {\widehat{\sigma }}_{fz}=\frac{s\widehat{E}\left(z,s\right)}{\left(1+\nu \right)\left(1-2\nu \right)}\left[\nu {\widehat{u}}_{f,x}+\nu {\widehat{v}}_{f,y}+\left(1-\nu \right){\widehat{w}}_{f,z}\right]\\ {\widehat{\tau }}_{\mathit{fxy}}=\frac{s\widehat{E}\left(z,s\right)}{2\left(1+\nu \right)}\left({\widehat{u}}_{f,y}+{\widehat{v}}_{f,x}\right)\\ {\widehat{\tau }}_{\mathit{fxz}}=\frac{s\widehat{E}\left(z,s\right)}{2\left(1+\nu \right)}\left({\widehat{u}}_{f,z}+{\widehat{w}}_{f,x}\right)\mathit{}\\ {\widehat{\tau }}_{\mathit{fyz}}=\frac{s\widehat{E}\left(z,s\right)}{2\left(1+\nu \right)}\left({\widehat{v}}_{f,z}+{\widehat{w}}_{f,y}\right),\end{array}$$(4) , state-space governing differential equations in matrix form are derived as: (5) $$\frac{\mathit{d}}{\mathit{\text{dz}}}{\mathit{\delta }}_{\mathit{f}}{=\mathit{G}}_{\mathit{f}}{\mathit{\delta }}_{\mathit{f}},$$(5) where ${\mathit{\delta }}_{\mathit{f}}={\left[{\widehat{\sigma }}_{\mathit{\text{fz}}}{,\widehat{u}}_{\mathit{f}}{,\widehat{v}}_{\mathit{f}}{,\widehat{w}}_{\mathit{f}}{,\widehat{\tau }}_{\mathit{\text{fxz}}}{,\widehat{\tau }}_{\mathit{\text{fyz}}}\right]}^{\mathit{T}}$ are the state variables and ${\mathit{G}}_{\mathit{f}}$ is the square coefficient matrix, defined in Eq. (A-2).

Additionally, induced variables (in-plane stresses) can be written in terms of the state variables as follow: (6) $$\begin{array}{c}{\widehat{\sigma }}_{\mathit{\text{fx}}}=\frac{1}{\left(1-\mathit{\nu }\right)}\left[{\nu \widehat{\mathit{\text{σ}}}}_{\mathit{\text{fz}}}+\frac{s\widehat{\mathit{\text{E}}}\left(\mathit{z},\mathit{s}\right)}{1+\mathit{\nu }}{\widehat{u}}_{\mathit{f},\mathit{x}}+\frac{s\widehat{\mathit{\text{E}}}\left(\mathit{z},\mathit{s}\right)}{1+\mathit{\nu }}{\widehat{v}}_{\mathit{f},\mathit{y}}\right],\\ {\widehat{\sigma }}_{\mathit{\text{fy}}}=\frac{1}{\left(1-\mathit{\nu }\right)}\left[{\nu \widehat{\mathit{\text{σ}}}}_{\mathit{\text{fz}}}+\frac{s\widehat{\mathit{\text{E}}}\left(\mathit{z},\mathit{s}\right)}{1+\mathit{\nu }}{\widehat{u}}_{\mathit{f},\mathit{x}}+\frac{s\widehat{\mathit{\text{E}}}\left(\mathit{z},\mathit{s}\right)}{1+\mathit{\nu }}{\widehat{v}}_{\mathit{f},\mathit{y}}\right],\\ {\widehat{\tau }}_{\mathit{\text{fxy}}}=\frac{s\widehat{\mathit{\text{E}}}\left(\mathit{z},\mathit{s}\right)}{2\left(1+\mathit{\nu }\right)}\left({\widehat{u}}_{\mathit{f},\mathit{y}}{+\widehat{\mathit{\text{u}}}}_{\mathit{f},\mathit{x}}\right).\end{array}$$(6)

In the time domain, the relaxation modulus, $\mathit{E},$ of the FGM layer is assumed in the form of a Prony series (Norouzi and Alibeigloo 2018), (7) $$\mathit{E}\left(\mathit{z},\mathit{t}\right){=\mathit{E}}_{1\mathit{m}}{\left(\frac{\mathit{z}}{{\mathit{h}}_{\mathit{p}}}\right)}^{{\mathit{m}}_1}{+\mathit{E}}_{2\mathit{m}}{\left(\frac{\mathit{z}}{{\mathit{h}}_{\mathit{p}}}\right)}^{{\mathit{m}}_2}{\mathit{e}}^{\left(-\frac{\mathit{t}}{\tau }\right)},$$(7) in which ${\mathit{m}}_1=\frac{\hspace{0.17em}\mathit{\text{ln}}\hspace{0.17em}\left(\frac{{\mathit{E}}_{1\mathit{c}}}{{\mathit{E}}_{1\mathit{m}}}\right)}{\hspace{0.17em}\mathit{\text{ln}}\hspace{0.17em}\left(\frac{{\mathit{h}}_{\mathit{f}}{+\mathit{h}}_{\mathit{p}}}{{\mathit{h}}_{\mathit{p}}}\right)}$ and ${\mathit{m}}_2=\frac{\hspace{0.17em}\mathit{\text{ln}}\hspace{0.17em}\left(\frac{{\mathit{E}}_{2\mathit{c}}}{{\mathit{E}}_{2\mathit{m}}}\right)}{\hspace{0.17em}\mathit{\text{ln}}\hspace{0.17em}\left(\frac{{\mathit{h}}_{\mathit{f}}{+\mathit{h}}_{\mathit{p}}}{{\mathit{h}}_{\mathit{p}}}\right)}.$

The subscripts $\mathit{m}$ and $\mathit{c}$ refer to metal and ceramic constituents in the FGM layer, which could be copper and tungsten or aluminum and alumina (Norouzi and Alibeigloo 2018). $\mathit{\tau }$ denotes the relaxation time constant, which is considered identical for the two constituents.

Piezoelectric layer

Both piezoelectric layers are assumed to exhibit linear piezoelastic material behavior; in particular, we only concern ourselves with orthorhombic piezoelectric materials of crystal class 2 mm, used for the actuator layers; and with tetragonal piezoelectric materials of crystal class 4 mm, for which the plane of the layer is the isotropic, used for the sensor layer. In these cases, the constitutive relations can be written as: (8) $$\left\{{\mathit{\sigma }}_{\mathit{p}}\right\}=\left[{\mathit{C}}_{\mathit{p}}\right]\left\{{\mathit{\epsilon }}_{\mathit{p}}\right\}{-\mathit{e}}^{\mathit{T}}{\mathit{E}}_{\mathit{p}}{,\mathit{}\mathit{D}=\mathit{e}\mathit{\epsilon }}_{\mathit{p}}{+\mathit{\eta }\mathit{E}}_{\mathit{p}},$$(8) where ${\mathit{C}}_{\mathit{p}}{,\mathit{\epsilon }}_{\mathit{p}}{,\mathit{\sigma }}_{\mathit{p}},$ ${\mathit{D},\mathit{E}}_{\mathit{p}},\mathit{\eta }$ and $\mathit{e}$ are defined in Eqs. (A-3)–(A-5).

The equilibrium equations, in the absence of body forces, are as (9) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{px}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{pxy}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{pxz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{pxy}},\mathit{x}}{+\mathit{\sigma }}_{\mathit{\text{py}},\mathit{y}}{+\mathit{\tau }}_{\mathit{\text{pyz}},\mathit{z}}=0\\ {\mathit{\tau }}_{\mathit{\text{pxz}},\mathit{x}}{+\mathit{\tau }}_{\mathit{\text{pyz}},\mathit{y}}{+\mathit{\sigma }}_{\mathit{\text{pz}},\mathit{z}}=0.\end{array}$$(9)

And, for the non-conducting piezoelectric layers, we also have the conservation of electrical charges (10) $${\mathit{D}}_{\mathit{x},\mathit{x}}{+\mathit{D}}_{\mathit{y},\mathit{y}}{+\mathit{D}}_{\mathit{z},\mathit{z}}=0.$$(10)

The relations between the electric field and electric potential are (11) $${\mathit{E}}_{\mathit{x}}{=-\mathit{\psi }}_{,\mathit{x}}{,\mathit{\text{E}}}_{\mathit{y}}{=-\mathit{\psi }}_{,\mathit{y}}{,\mathit{\text{E}}}_{\mathit{z}}{=-\mathit{\psi }}_{,\mathit{z}}.$$(11)

At the upper surface of the actuator and the lower surface of the sensor, there are the following mechanical and electrical surface boundary conditions: (12) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{pz}}}=-\mathit{q},\hspace{1em}{\mathit{\tau }}_{\mathit{\text{pzx}}}{=\mathit{\tau }}_{\mathit{\text{pzy}}}=0,\hspace{1em}{\mathit{\psi }=\mathit{V}}_0,\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{z}=\mathit{h}\\ {\mathit{\sigma }}_{\mathit{\text{pz}}}{=\mathit{\tau }}_{\mathit{\text{pzx}}}{=\mathit{\tau }}_{\mathit{\text{pzy}}}=0,\hspace{1em}{\mathit{D}}_{\mathit{z}}=0,\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{z}=0.\end{array}$$(12)

Moreover, at the interfaces of the actuator and sensor with the viscoelastic layers, electric potential, $\mathit{\psi },$ is also zero.

At the boundaries, for a simply supported plate, the following conditions are satisfied: (13) $$\begin{array}{c}{\mathit{w}}_{\mathit{p}}{=\mathit{v}}_{\mathit{p}}{=0,\mathit{\text{σ}}}_{\mathit{\text{px}}}=0,\mathit{}\mathit{\psi }=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{x}=0,\mathit{a},\\ {\mathit{w}}_{\mathit{p}}{=\mathit{u}}_{\mathit{p}}{=0,\mathit{\text{σ}}}_{\mathit{\text{py}}}=0,\mathit{}\mathit{\psi }=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{y}=0,\mathit{b}.\end{array}$$(13)

Combining Eqs. (8)–(13), similar to the FGM layer, constitutive equations for the piezoelectric layers in terms of displacement components are derived; and with the similar manipulation, the following state-space governing equation in Laplace domain for the piezoelectric layers is yielded. (14) $$\frac{\mathit{d}}{\mathit{\text{dz}}}{\mathit{\delta }}_{\mathit{P}}{=\mathit{G}}_{\mathit{P}}{\mathit{\delta }}_{\mathit{P}},$$(14) where ${\delta }_P={\left\{\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cc}{\widehat{\sigma }}_{pz}& {\widehat{u}}_p\end{array}& {\widehat{v}}_p& {\widehat{w}}_p& {\widehat{\tau }}_{\mathit{pzx}}\end{array}& {\widehat{\tau }}_{\mathit{pzy}}& {\widehat{D}}_z& \widehat{\psi }\end{array}\right\}}^T;$ ${\mathit{G}}_{\mathit{P}}$ is a constant coefficient square matrix, defined in Eq. (A-6)(A-6) $$G_P=\left[\begin{array}{cccccccc}0& 0& 0& 0& -\frac{\partial }{\partial x}& -\frac{\partial }{\partial y}& 0& 0\\ 0& 0& 0& -\frac{\partial }{\partial x}& \frac{1}{Q_{55}}& 0& 0& -\frac{e_5}{Q_{55}}\frac{\partial }{\partial x}\\ 0& 0& 0& -\frac{\partial }{\partial y}& 0& \frac{1}{Q_{44}}& 0& -\frac{e_4}{Q_{44}}\frac{\partial }{\partial y}\\ \frac{{\eta }_3}{Q_{33}{\eta }_3+e_3^2}& -\frac{Q_{13}{\eta }_3+e_1{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial x}& G_{P43}\frac{\partial }{\partial y}& 0& 0& 0& G_{P47}& 0\\ -\frac{Q_{13}{\eta }_3+e_1{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial x}& G_{P52}& G_{P53}& 0& 0& 0& \frac{Q_{33}{e}_1-C_{13}{e}_3}{Q_{33}{\eta }_3+e_3^2}& 0\\ -\frac{Q_{23}{\eta }_3+e_2{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial y}& G_{P62}& G_{P63}& 0& 0& 0& \frac{Q_{33}{e}_1-Q_{23}{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial y}& 0\\ 0& 0& 0& 0& -\frac{e_5}{Q_{55}}\frac{\partial }{\partial x}& -\frac{e_4}{Q_{44}}\frac{\partial }{\partial y}& 0& \left(\frac{e_5^2}{Q_{55}}+{\eta }_1\right)\frac{{\partial }^2}{\partial x^2}+\left(\frac{e_4^2}{Q_{44}}+{\eta }_2\right)\frac{{\partial }^2}{\partial y^2}\\ \frac{e_3}{Q_{33}{\eta }_3+e_3^2}& -\left(\frac{Q_{13}{e}_3-Q_{33}{e}_1}{Q_{33}{\eta }_3+e_3^2}\right)\frac{\partial }{\partial x}& -\left(\frac{Q_{23}{e}_3-Q_{33}{e}_2}{Q_{33}{\eta }_3+e_3^2}\right)\frac{\partial }{\partial y}& 0& 0& 0& -\frac{Q_{33}}{Q_{33}{\eta }_3+e_3^2}& 0\end{array}\right]$$(A-6) .

Solution procedure

FGM layer

Relations for simply supported boundary conditions of the rectangular plate are (15) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{fx}}}{=0,\mathit{}\mathit{v}}_f{=\mathit{w}}_f=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{x}=0,\mathit{a},\\ {\mathit{\sigma }}_{\mathit{\text{fy}}}{=0,\mathit{}\mathit{u}}_f{=\mathit{w}}_f=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{y}=0,\mathit{b}.\end{array}$$(15)

The following Fourier series expansions for stress and the displacement components satisfy simply supported boundary conditions relation, Eq. (15)(15) $$\begin{array}{c}{\mathit{\sigma }}_{\mathit{\text{fx}}}{=0,\mathit{}\mathit{v}}_f{=\mathit{w}}_f=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{x}=0,\mathit{a},\\ {\mathit{\sigma }}_{\mathit{\text{fy}}}{=0,\mathit{}\mathit{u}}_f{=\mathit{w}}_f=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{y}=0,\mathit{b}.\end{array}$$(15) (16) $$\begin{array}{c}{\widehat{u}}_{\mathit{i}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{u}}_{\mathit{i}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{v}}_{\mathit{i}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{v}}_{\mathit{i}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\mathit{}\mathit{}\\ {\widehat{w}}_{\mathit{i}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{w}}_{\mathit{i}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\sigma }}_{\mathit{\text{iz}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\sigma }}_{\mathit{\text{iz}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\mathit{}\mathit{}\\ {\widehat{\tau }}_{\mathit{\text{iyz}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\tau }}_{\mathit{\text{iyz}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\tau }}_{\mathit{\text{ixz}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\tau }}_{\mathit{\text{ixz}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\sigma }}_{\mathit{\text{ix}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\sigma }}_{\mathit{\text{ix}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\sigma }}_{\mathit{\text{iy}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\sigma }}_{\mathit{\text{iy}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\tau }}_{\mathit{\text{ixy}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\tau }}_{\mathit{\text{ixy}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\mathit{}\mathit{}\\ \mathit{i}=\mathit{f},\mathit{p}\end{array}$$(16)

For convenience, the following non-dimensional quantities are introduced: (17) $$\begin{array}{l}({\overline{u}}_{\mathit{i}}{,\overline{v}}_{\mathit{i}}{,\overline{w}}_{\mathit{i}})=({\widehat{u}}_{\mathit{i}}^{\prime }{,\widehat{v}}_{\mathit{i}}^{\prime }{,\widehat{w}}_{\mathit{i}}^{\prime })\frac{1}{{(\widehat{q}/\mathit{\text{k}}}^{*})\mathit{h}},\mathit{}\overline{E}=\frac{\widehat{E}}{\widehat{q}},\\ ({\overline{\sigma }}_{\mathit{\text{ix}}}{,\overline{\sigma }}_{\mathit{\text{iy}}}{,\overline{\sigma }}_{\mathit{\text{iz}}}{,\overline{\tau }}_{\mathit{\text{izx}}}{,\overline{\tau }}_{\mathit{\text{izx}}}{,\overline{\tau }}_{\mathit{\text{ixy}}})=({\widehat{\sigma }}_{\mathit{\text{ix}}}^{\prime }{,\widehat{\sigma }}_{\mathit{\text{iy}}}^{\prime }{,\widehat{\sigma }}_{\mathit{\text{iz}}}^{\prime }{,\widehat{\tau }}_{\mathit{\text{izx}}}^{\prime }{,\widehat{\tau }}_{\mathit{\text{izy}}}^{\prime }{,\widehat{\tau }}_{\mathit{\text{ixy}}}^{\prime })\frac{1}{\widehat{q}},\\ \overline{\mathit{\text{x}}}=\frac{\mathit{x}}{\mathit{h}}\mathit{},\mathit{}\overline{y}=\frac{\mathit{y}}{\mathit{h}}\mathit{},\mathit{}\overline{z}=\frac{\mathit{z}}{\mathit{h}}\mathit{},\\ {\overline{P}}_{\mathit{m}}{=\mathit{\text{aP}}}_{\mathit{m}}{,\mathit{}\overline{P}}_{\mathit{n}}{=\mathit{\text{bP}}}_{\mathit{n}},\\ \mathit{i}=\mathit{f},\mathit{p}.\end{array}$$(17)

Here ${\mathit{k}}^{*}=1\mathit{\text{GPa}}$ is the scale factor and $\widehat{q}$ is the applied mechanical load.

Then, we derive a non-dimensional version of the state-space formulation introduced in Eq. (5)(5) $$\frac{\mathit{d}}{\mathit{\text{dz}}}{\mathit{\delta }}_{\mathit{f}}{=\mathit{G}}_{\mathit{f}}{\mathit{\delta }}_{\mathit{f}},$$(5) for the non-dimensional expansion coefficients of the Fourier series as (18) $$\frac{\mathit{d}}{d\overline{\mathit{\text{z}}}}{\overline{\delta }}_{\mathit{f}}{=\overline{G}}_{\mathit{f}}{\overline{\delta }}_{\mathit{f}},$$(18) where ${\overline{\delta }}_f={\left[{\overline{\sigma }}_{fz},{\overline{u}}_f,{\overline{v}}_f,{\overline{w}}_f,{\overline{\tau }}_{\mathit{fxz}},{\overline{\tau }}_{\mathit{fyz}}\right]}^T;$ ${\overline{G}}_{\mathit{f}}$ is the non-dimensional square matrix of coefficients, defined in Eq. (A-7)(A-7) $${\overline{G}}_f=\left[\begin{array}{cccccc}0& 0& 0& 0& \left(\frac{h}{a}\right){\overline{P}}_m& \left(\frac{h}{b}\right){\overline{P}}_n\\ 0& 0& 0& -\left(\frac{h}{a}\right){\overline{P}}_m& \frac{2\left(1+\nu \right)}{s\overline{E}}& 0\\ 0& 0& 0& -\left(\frac{h}{b}\right){\overline{P}}_n& 0& \frac{2\left(1+\nu \right)}{s\overline{E}}\\ \frac{\left(1+\nu \right)\left(1-2\nu \right)}{\left(1-\nu \right)s\overline{E}}& {\overline{P}}_m\left(\frac{h}{a}\right)\frac{\nu }{\left(1-\nu \right)}& {\overline{P}}_n\left(\frac{h}{b}\right)\frac{\nu }{\left(1-\nu \right)}& 0& 0& 0\\ -{\overline{P}}_m\left(\frac{h}{a}\right)\frac{\nu }{\left(1-\nu \right)}& {\overline{P}}^2{}_m{\left(\frac{h}{a}\right)}^2\frac{s\overline{E}}{\left(1-{\nu }^2\right)}-{\overline{P}}^2{}_n{\left(\frac{h}{b}\right)}^2\frac{s\overline{E}}{2\left(1+\nu \right)}& {\overline{P}}_m{\overline{P}}_n\left(\frac{h}{a}\right)\left(\frac{h}{b}\right)\frac{s\overline{E}}{2\left(1+\nu \right)}& 0& 0& 0\\ -{\overline{P}}_n\left(\frac{h}{b}\right)\frac{\nu }{\left(1-\nu \right)}& {\overline{P}}_m{\overline{P}}_n\left(\frac{h}{a}\right)\left(\frac{h}{b}\right)\frac{s\overline{E}}{2\left(1+\nu \right)}& {\overline{P}}^2{}_n{\left(\frac{h}{b}\right)}^2\frac{s\overline{E}}{\left(1-{\nu }^2\right)}+{\overline{P}}^2{}_m{\left(\frac{h}{a}\right)}^2\frac{s\overline{E}}{2\left(1+\nu \right)}& 0& 0& 0\end{array}\right]$$(A-7) .

Also, the non-dimensional state variables can be used to derive the non-dimensional induced variables: (19) $$\left\{\begin{array}{c}{\overline{\sigma }}_{\mathit{\text{fx}}}\\ {\overline{\sigma }}_{\mathit{\text{fy}}}\\ {\overline{\tau }}_{\mathit{\text{xy}}}\end{array}\right\}=\left[\overline{F}\right]\left\{\begin{array}{c}{\overline{\sigma }}_{\mathit{\text{fz}}}\\ {\overline{u}}_{\mathit{f}}\\ {\overline{v}}_{\mathit{f}}\end{array}\right\},$$(19) where $\overline{F}$ is a coefficient matrix, defined in Eq. (A-8)(8) $$\left\{{\mathit{\sigma }}_{\mathit{p}}\right\}=\left[{\mathit{C}}_{\mathit{p}}\right]\left\{{\mathit{\epsilon }}_{\mathit{p}}\right\}{-\mathit{e}}^{\mathit{T}}{\mathit{E}}_{\mathit{p}}{,\mathit{}\mathit{D}=\mathit{e}\mathit{\epsilon }}_{\mathit{p}}{+\mathit{\eta }\mathit{E}}_{\mathit{p}},$$(8) .

Since the transfer matric ${\overline{G}}_{\mathit{f}}$ is not constant because of the variation of the material properties along the thickness direction, we employ a layer-wise approximation dividing the FGM layer into ${\mathit{N}}_{\mathit{f}}$ fictitious thin layers whose material parameters are assumed constant along the z-direction; hence, the full transfer matrix for the fictitious layers ${\overline{G}}_{\mathit{k}},$ can be considered constant in each fictitious layer, $\mathit{k}=1,\dots {,\mathit{N}}_{\mathit{f}},$ leading to the following solution to Eq. (18)(18) $$\frac{\mathit{d}}{d\overline{\mathit{\text{z}}}}{\overline{\delta }}_{\mathit{f}}{=\overline{G}}_{\mathit{f}}{\overline{\delta }}_{\mathit{f}},$$(18) (20) $$\overline{\delta }\left(\overline{z}\right){=\mathit{e}}^{{\overline{G}}_{\mathit{k}}\overline{z}}{\mathit{e}}^{{-\overline{G}}_{\mathit{k}}{\overline{z}}_{\mathit{k}-1}}{\overline{\delta }}_{0\mathit{k}},\hspace{1em}{\overline{z}}_{\mathit{k}-1}\le \overline{z}\le {\overline{z}}_{\mathit{k}},$$(20) in which ${\overline{z}}_{\mathit{k}-1}{=\overline{h}}_{\mathit{p}}+\frac{(\mathit{k}-1){\overline{h}}_{\mathit{f}}}{{\mathit{N}}_{\mathit{f}}},{\overline{z}}_{\mathit{k}}{=\overline{h}}_{\mathit{p}}+\frac{{k\overline{\mathit{\text{h}}}}_{\mathit{f}}}{{\mathit{N}}_{\mathit{f}}},{\mathit{\delta }}_{0\mathit{k}}{=\mathit{\delta }}_{\mathit{k}}({\overline{z}=\overline{z}}_{\mathit{k}-1}).\hspace{1em}$

Then Eq. (13)(13) $$\begin{array}{c}{\mathit{w}}_{\mathit{p}}{=\mathit{v}}_{\mathit{p}}{=0,\mathit{\text{σ}}}_{\mathit{\text{px}}}=0,\mathit{}\mathit{\psi }=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{x}=0,\mathit{a},\\ {\mathit{w}}_{\mathit{p}}{=\mathit{u}}_{\mathit{p}}{=0,\mathit{\text{σ}}}_{\mathit{\text{py}}}=0,\mathit{}\mathit{\psi }=0\hspace{1em}\mathit{\text{at}}\hspace{1em}\mathit{y}=0,\mathit{b}.\end{array}$$(13) at ${\overline{z}}_{\mathit{k}}$ becomes (21) $${\overline{\delta }}_{\mathit{k}}\left({\overline{z}}_{\mathit{k}}\right){=\mathit{M}}_{\mathit{k}}{\overline{\delta }}_{0\mathit{k}},$$(21) where ${\mathit{M}}_{\mathit{k}}{=\mathit{e}}^{({\overline{G}}_{\mathit{k}}\frac{{\overline{h}}_{\mathit{f}}}{{\mathit{N}}_{\mathit{f}}})}.$

Regarding the continuity condition of the state variables at the interface of each two adjacent fictitious layers, state variables at the top and the bottom surfaces of the viscoelastic layer are related as follows: (22) $${\overline{\delta }}_{\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right)}{=\mathit{M}}_{\mathit{f}}{\overline{\delta }}_0,$$(22) where ${\mathit{M}}_{\mathit{f}}={\prod }_{\mathit{k}=\mathit{N}}^1{\mathit{e}}^{({\overline{G}}_{\mathit{\text{fk}}}\frac{{\overline{h}}_{\mathit{f}}}{{\mathit{N}}_{\mathit{f}}})};$ ${\overline{\delta }}_0$ and ${\overline{\delta }}_{((h_p+h_f)/h)}$ are the state vector at the bottom and top surfaces of the viscoelastic layer, respectively.

Piezoelectric layer

The following dimensionless quantities are defined for the piezoelectric layers: (23) $$\begin{array}{l}{(\overline{E}}_{\mathit{x}}{,\overline{E}}_{\mathit{y}}{,\overline{E}}_{\mathit{z}}{)=(\mathit{E}}_{\mathit{x}}{,\mathit{E}}_{\mathit{y}}{,\mathit{E}}_{\mathit{z}})\frac{|{\mathit{d}}_1|}{{(\widehat{q}/\mathit{\text{k}}}^{*})},({\overline{D}}_{\mathit{x}}{,\overline{D}}_{\mathit{y}}{,\overline{D}}_{\mathit{z}})=({\widehat{D}}^{\prime }{}_{\mathit{x}}{,\widehat{D}}^{\prime }{}_{\mathit{y}}{,\widehat{D}}^{\prime }{}_{\mathit{z}})\frac{1}{\mathit{q}|{\mathit{d}}_1|},\\ ({\overline{\eta }}_1{,\overline{\eta }}_2{,\overline{\eta }}_3)=({\mathit{\eta }}_1{,\mathit{\eta }}_2{,\mathit{\eta }}_3)\frac{1}{{\mathit{k}}^{*}{|{\mathit{d}}_1|}^2},\hspace{1em}({\overline{e}}_1{,\overline{e}}_2{,\overline{e}}_3)=({\mathit{e}}_1{,\mathit{e}}_2{,\mathit{e}}_3)\frac{1}{{\mathit{k}}^{*}|{\mathit{d}}_1|},\\ {\widehat{\psi }}^{\prime }=\frac{|{\mathit{d}}_1|\overline{\psi }}{{(\widehat{q}/\mathit{\text{k}}}^{*})\mathit{h}},\hspace{1em}{\overline{C}}_{\mathit{\text{ij}}}=\frac{{\mathit{C}}_{\mathit{\text{ij}}}}{\mathit{k}*}.\end{array}$$(23)

In addition to the Fourier series solution used in Eq. (17)(17) $$\begin{array}{l}({\overline{u}}_{\mathit{i}}{,\overline{v}}_{\mathit{i}}{,\overline{w}}_{\mathit{i}})=({\widehat{u}}_{\mathit{i}}^{\prime }{,\widehat{v}}_{\mathit{i}}^{\prime }{,\widehat{w}}_{\mathit{i}}^{\prime })\frac{1}{{(\widehat{q}/\mathit{\text{k}}}^{*})\mathit{h}},\mathit{}\overline{E}=\frac{\widehat{E}}{\widehat{q}},\\ ({\overline{\sigma }}_{\mathit{\text{ix}}}{,\overline{\sigma }}_{\mathit{\text{iy}}}{,\overline{\sigma }}_{\mathit{\text{iz}}}{,\overline{\tau }}_{\mathit{\text{izx}}}{,\overline{\tau }}_{\mathit{\text{izx}}}{,\overline{\tau }}_{\mathit{\text{ixy}}})=({\widehat{\sigma }}_{\mathit{\text{ix}}}^{\prime }{,\widehat{\sigma }}_{\mathit{\text{iy}}}^{\prime }{,\widehat{\sigma }}_{\mathit{\text{iz}}}^{\prime }{,\widehat{\tau }}_{\mathit{\text{izx}}}^{\prime }{,\widehat{\tau }}_{\mathit{\text{izy}}}^{\prime }{,\widehat{\tau }}_{\mathit{\text{ixy}}}^{\prime })\frac{1}{\widehat{q}},\\ \overline{\mathit{\text{x}}}=\frac{\mathit{x}}{\mathit{h}}\mathit{},\mathit{}\overline{y}=\frac{\mathit{y}}{\mathit{h}}\mathit{},\mathit{}\overline{z}=\frac{\mathit{z}}{\mathit{h}}\mathit{},\\ {\overline{P}}_{\mathit{m}}{=\mathit{\text{aP}}}_{\mathit{m}}{,\mathit{}\overline{P}}_{\mathit{n}}{=\mathit{\text{bP}}}_{\mathit{n}},\\ \mathit{i}=\mathit{f},\mathit{p}.\end{array}$$(17) , the following series expansion are used for simply supported boundary conditions in the piezoelectric layer: (24) $$\begin{array}{c}\widehat{\psi }=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\psi }}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right),\\ {\widehat{D}}_{\mathit{z}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{D}}^{\prime }{}_{\mathit{z}}\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right).\end{array}$$(24)

By substituting Eqs. (16)(16) $$\begin{array}{c}{\widehat{u}}_{\mathit{i}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{u}}_{\mathit{i}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{v}}_{\mathit{i}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{v}}_{\mathit{i}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\mathit{}\mathit{}\\ {\widehat{w}}_{\mathit{i}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{w}}_{\mathit{i}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\sigma }}_{\mathit{\text{iz}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\sigma }}_{\mathit{\text{iz}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\mathit{}\mathit{}\\ {\widehat{\tau }}_{\mathit{\text{iyz}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\tau }}_{\mathit{\text{iyz}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\tau }}_{\mathit{\text{ixz}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\tau }}_{\mathit{\text{ixz}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\sigma }}_{\mathit{\text{ix}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\sigma }}_{\mathit{\text{ix}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\sigma }}_{\mathit{\text{iy}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\sigma }}_{\mathit{\text{iy}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\\ {\widehat{\tau }}_{\mathit{\text{ixy}}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\tau }}_{\mathit{\text{ixy}}}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{cos}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right)\mathit{}\mathit{}\\ \mathit{i}=\mathit{f},\mathit{p}\end{array}$$(16) , (23)(23) $$\begin{array}{l}{(\overline{E}}_{\mathit{x}}{,\overline{E}}_{\mathit{y}}{,\overline{E}}_{\mathit{z}}{)=(\mathit{E}}_{\mathit{x}}{,\mathit{E}}_{\mathit{y}}{,\mathit{E}}_{\mathit{z}})\frac{|{\mathit{d}}_1|}{{(\widehat{q}/\mathit{\text{k}}}^{*})},({\overline{D}}_{\mathit{x}}{,\overline{D}}_{\mathit{y}}{,\overline{D}}_{\mathit{z}})=({\widehat{D}}^{\prime }{}_{\mathit{x}}{,\widehat{D}}^{\prime }{}_{\mathit{y}}{,\widehat{D}}^{\prime }{}_{\mathit{z}})\frac{1}{\mathit{q}|{\mathit{d}}_1|},\\ ({\overline{\eta }}_1{,\overline{\eta }}_2{,\overline{\eta }}_3)=({\mathit{\eta }}_1{,\mathit{\eta }}_2{,\mathit{\eta }}_3)\frac{1}{{\mathit{k}}^{*}{|{\mathit{d}}_1|}^2},\hspace{1em}({\overline{e}}_1{,\overline{e}}_2{,\overline{e}}_3)=({\mathit{e}}_1{,\mathit{e}}_2{,\mathit{e}}_3)\frac{1}{{\mathit{k}}^{*}|{\mathit{d}}_1|},\\ {\widehat{\psi }}^{\prime }=\frac{|{\mathit{d}}_1|\overline{\psi }}{{(\widehat{q}/\mathit{\text{k}}}^{*})\mathit{h}},\hspace{1em}{\overline{C}}_{\mathit{\text{ij}}}=\frac{{\mathit{C}}_{\mathit{\text{ij}}}}{\mathit{k}*}.\end{array}$$(23) , and (24)(24) $$\begin{array}{c}\widehat{\psi }=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{\psi }}^{\prime }\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right),\\ {\widehat{D}}_{\mathit{z}}=\sum _{\mathit{m}=0}^{\infty }\sum _{\mathit{n}=0}^{\infty }{\widehat{D}}^{\prime }{}_{\mathit{z}}\mathit{\text{(z)}}\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{m}}\mathit{x}\right)\hspace{0.17em}\mathit{\text{sin}}\hspace{0.17em}\left({\mathit{p}}_{\mathit{n}}\mathit{y}\right).\end{array}$$(24) into Eq. (18)(18) $$\frac{\mathit{d}}{d\overline{\mathit{\text{z}}}}{\overline{\delta }}_{\mathit{f}}{=\overline{G}}_{\mathit{f}}{\overline{\delta }}_{\mathit{f}},$$(18) , the following non-dimensional state-space equation in matrix form for the piezoelectric layers is derived: (25-a) $$\frac{\mathit{d}}{\mathit{\text{dz}}}{\overline{\delta }}_{\mathit{P}}{=\overline{G}}_{\mathit{P}}{\overline{\delta }}_{\mathit{P}},$$(25-a) where ${\overline{\delta }}_P={\left\{\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cc}{\overline{\sigma }}_{pz}& {\overline{u}}_p\end{array}& {\overline{v}}_p{}_{}& {\overline{w}}_p& {\overline{\tau }}_{\mathit{pzx}}\end{array}& {\overline{\tau }}_{\mathit{pzy}}& {\overline{D}}_z& \overline{\psi }\end{array}\right\}}^T;$ the matrix ${\overline{G}}_{\mathit{P}}$ is defined in Eq. (A-9). In term of state variables, the induced variables for the piezoelectric layers can be derived as follow: (25-b) $$\left\{\begin{array}{c}{\overline{\sigma }}_{\mathit{\text{px}}}\\ {\overline{\sigma }}_{\mathit{\text{py}}}\\ {\overline{\tau }}_{\mathit{\text{pxy}}}\end{array}\right\}=\left[{\overline{F}}_{\mathit{p}}\right]\left\{\begin{array}{c}{\overline{\sigma }}_{\mathit{\text{pz}}}\\ {\overline{u}}_{\mathit{p}}\\ {\overline{v}}_{\mathit{p}}\\ {\overline{D}}_{\mathit{z}}\end{array}\right\},$$(25-b) where ${\overline{F}}_{\mathit{p}}$ is the coefficient matrix, defined in Eq. (A-10).

In the sensor and actuator layers with simply supported boundary conditions, the general solution for Eq. (25(25-a) $$\frac{\mathit{d}}{\mathit{\text{dz}}}{\overline{\delta }}_{\mathit{P}}{=\overline{G}}_{\mathit{P}}{\overline{\delta }}_{\mathit{P}},$$(25-a) -a), respectively, are (26-a) $${\overline{\delta }}_{\mathit{s}}\left(\overline{z}\right){=\mathit{e}}^{{\overline{G}}_{\mathit{z}}\overline{z}}\overline{\delta }\left(0\right)\hspace{1em}at\hspace{1em}0\le \overline{z}\le {\mathit{h}}_{\mathit{p}}/\mathit{h}$$(26-a) (26-b) $${\overline{\delta }}_{\mathit{a}}\left(\overline{z}\right){=\mathit{e}}^{{\overline{G}}_{\mathit{a}}\overline{z}}{\mathit{e}}^{{-\overline{G}}_{\mathit{a}}{(\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}})/\mathit{h}}\overline{\delta }\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right)\hspace{1em}\mathit{\text{at}}\hspace{1em}\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\le \overline{z}\le \mathit{h}/\mathit{h}.$$(26-b)

Equations (26(26-a) $${\overline{\delta }}_{\mathit{s}}\left(\overline{z}\right){=\mathit{e}}^{{\overline{G}}_{\mathit{z}}\overline{z}}\overline{\delta }\left(0\right)\hspace{1em}at\hspace{1em}0\le \overline{z}\le {\mathit{h}}_{\mathit{p}}/\mathit{h}$$(26-a) -a) and (26-b) at the above surfaces of actuator and sensor layers are, respectively, (27-a) $${\overline{\delta }}_{\mathit{a}}\left(\mathit{h}/\mathit{h}\right){=\mathit{M}}^{\mathit{a}}{\overline{\delta }}_{\mathit{a}}\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right).$$(27-a) (27-b) $${\overline{\delta }}_{\mathit{s}}\left({\mathit{h}}_{\mathit{p}}\right){=\mathit{M}}^{\mathit{s}}{\overline{\delta }}_{\mathit{s}}\left(0\right),$$(27-b) where ${\mathit{M}}^{\mathit{a}}{=\mathit{e}}^{\left[{\overline{G}}_{\mathit{a}}{\mathit{h}}_{\mathit{p}}/\mathit{h}\right]},$${\mathit{M}}^{\mathit{s}}{=\mathit{e}}^{\left[{\overline{G}}_{\mathit{s}}{\mathit{h}}_{\mathit{p}}/\mathit{h}\right]}.$

Based on Eq. (27a), at the lower surface of the actuator layer, electric displacement is (Feri, Alibeigloo, and Pasha Zanoosi, 2016): (28) $$\begin{array}{c}{\overline{D}}_{\mathit{z}}\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right)=\frac{1}{{\mathit{m}}_{87}^{\mathit{a}}}\left[\overline{\psi }\left(\mathit{h}/\mathit{h}\right)-\left[{\mathit{m}}_{8\mathit{j}}^{\mathit{a}}\right]{\overline{\delta }}^{\mathit{m}}\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right)\right],\\ \mathit{j}=1,\dots ,6.\end{array}$$(28) where ${\mathit{m}}_{8\mathit{j}}^{\mathit{a}}$ are the $8\mathrm{th}$ row and $\mathit{j}\mathrm{th}$ column elements of matrix $M^a$ and ${\overline{\delta }}^m((h_p+h_f)/h)$ is the mechanical part of the state variables of the actuator at its interface with the viscoelastic layer, ${\overline{\delta }}_p{}^m={\left\{{\overline{\sigma }}_{pz},{\overline{u}}_p,{\overline{v}}_p,{\overline{w}}_p,{\overline{\tau }}_{\mathit{pzx}},{\overline{\tau }}_{\mathit{pzy}}\right\}}^T$ inserting Eq. (27(27-a) $${\overline{\delta }}_{\mathit{a}}\left(\mathit{h}/\mathit{h}\right){=\mathit{M}}^{\mathit{a}}{\overline{\delta }}_{\mathit{a}}\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right).$$(27-a) -b) and relation $\overline{\psi }((h_p+h_f)/h)=0$ into Eq. (26(26-a) $${\overline{\delta }}_{\mathit{s}}\left(\overline{z}\right){=\mathit{e}}^{{\overline{G}}_{\mathit{z}}\overline{z}}\overline{\delta }\left(0\right)\hspace{1em}at\hspace{1em}0\le \overline{z}\le {\mathit{h}}_{\mathit{p}}/\mathit{h}$$(26-a) -a) yields: (29) $${\overline{\delta }}^{\mathit{m}}\left(\mathit{h}/\mathit{h}\right){=\mathit{T}}_{\mathit{a}}{\overline{\delta }}^{\mathit{m}}\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right){+\mathit{C}}_{\mathit{a}},$$(29) where ${\mathit{T}}_{\mathit{a}}=\left[\left[{\mathit{m}}_{\mathit{\text{ij}}}^{\mathit{a}}\right]-\left\{{\mathit{m}}_{\mathit{i}7}^{\mathit{a}}\right\}\frac{1}{{\mathit{m}}_{87}^{\mathit{a}}}\left[{\mathit{m}}_{8\mathit{j}}^{\mathit{a}}\right]\right],\hspace{1em}{\mathit{C}}_{\mathit{a}}=\left\{{\mathit{m}}_{\mathit{i}7}^{\mathit{a}}\right\}\frac{\overline{\psi }(\mathit{h})}{{\mathit{m}}_{87}^{\mathit{a}}},\hspace{1em}\mathit{i}=1,2,\dots ,6,\mathit{}\mathit{j}=1,2,\dots ,6.$

With the aid of Eq. (27b), at the lower surface of the sensor, the electric potential, $\overline{\psi },$ can be derived as: (30) $$\overline{\psi }\left(0\right)=\frac{1}{{\mathit{m}}_{88}^{\mathit{s}}}\left(\left[{\mathit{m}}_{8\mathit{j}}^{\mathit{s}}\right]{\overline{\delta }}^{\mathit{m}}\left(0\right)\right),\mathit{}\mathit{j}=1,2,\dots ,6.$$(30)

Setting ${\mathit{D}}_{\mathit{z}}=0$ and substituting Eq. (30)(30) $$\overline{\psi }\left(0\right)=\frac{1}{{\mathit{m}}_{88}^{\mathit{s}}}\left(\left[{\mathit{m}}_{8\mathit{j}}^{\mathit{s}}\right]{\overline{\delta }}^{\mathit{m}}\left(0\right)\right),\mathit{}\mathit{j}=1,2,\dots ,6.$$(30) into Eq. (27(27-a) $${\overline{\delta }}_{\mathit{a}}\left(\mathit{h}/\mathit{h}\right){=\mathit{M}}^{\mathit{a}}{\overline{\delta }}_{\mathit{a}}\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right).$$(27-a) -b) results in (31) $${\overline{\delta }}_{\mathit{p}}^{\mathit{m}}\left({\mathit{h}}_{\mathit{p}}/\mathit{h}\right){=\mathit{T}}_{\mathit{s}}{\overline{\delta }}_{\mathit{p}}^{\mathit{m}}\left(0\right),$$(31) where $${\mathit{T}}_{\mathit{s}}=\left[\left[{\mathit{m}}_{\mathit{i},\mathit{j}}^{\mathit{s}}\right]-\left\{{\mathit{m}}_{\mathit{i},8}^{\mathit{s}}\right\}\frac{1}{{\mathit{m}}_{8,8}^{\mathit{s}}}\left[{\mathit{m}}_{8,\mathit{j}}^{\mathit{s}}\right]\right],\hspace{1em}\mathit{i}=1,2,\dots ,6\hspace{1em}\mathit{j}=1,2,\dots ,6.$$

Because the mechanical tractions and displacements at the interfaces of the actuator, FGM, and sensor are continuous, using Eqs. (22)(22) $${\overline{\delta }}_{\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right)}{=\mathit{M}}_{\mathit{f}}{\overline{\delta }}_0,$$(22) , (29)(29) $${\overline{\delta }}^{\mathit{m}}\left(\mathit{h}/\mathit{h}\right){=\mathit{T}}_{\mathit{a}}{\overline{\delta }}^{\mathit{m}}\left(\left({\mathit{h}}_{\mathit{p}}{+\mathit{h}}_{\mathit{f}}\right)/\mathit{h}\right){+\mathit{C}}_{\mathit{a}},$$(29) , and (31)(31) $${\overline{\delta }}_{\mathit{p}}^{\mathit{m}}\left({\mathit{h}}_{\mathit{p}}/\mathit{h}\right){=\mathit{T}}_{\mathit{s}}{\overline{\delta }}_{\mathit{p}}^{\mathit{m}}\left(0\right),$$(31) , the relationship between the mechanical state variables at the top and bottom surfaces of the plate follows as: (32) $${\overline{\delta }}^{\mathit{m}}\left(\mathit{h}\right){=\mathit{S}\overline{\delta }}^{\mathit{m}}\left(0\right){+\mathit{C}}_{\mathit{a}},$$(32) where ${\mathit{S}=\mathit{T}}_{\mathit{a}}{\mathit{\text{MT}}}_{\mathit{s}}.$

Equation (32)(32) $${\overline{\delta }}^{\mathit{m}}\left(\mathit{h}\right){=\mathit{S}\overline{\delta }}^{\mathit{m}}\left(0\right){+\mathit{C}}_{\mathit{a}},$$(32) reduces to the following form by using surface traction at the bottom and top surfaces of the plate. (33) $$\left[\begin{array}{ccc}{\mathit{s}}_{12}& {\mathit{s}}_{13}& {\mathit{s}}_{14}\\ {\mathit{s}}_{52}& {\mathit{s}}_{53}& {\mathit{s}}_{54}\\ {\mathit{s}}_{62}& {\mathit{s}}_{63}& {\mathit{s}}_{64}\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_0\\ {\overline{v}}_0\\ {\overline{w}}_0\end{array}\right\}=\left\{\begin{array}{c}{\mathit{C}}_{\mathit{a}}\left(1\right)\\ {\mathit{C}}_{\mathit{a}}\left(5\right)\\ {\mathit{C}}_{\mathit{a}}\left(6\right)\end{array}\right\}.$$(33)

From Eq. (33)(33) $$\left[\begin{array}{ccc}{\mathit{s}}_{12}& {\mathit{s}}_{13}& {\mathit{s}}_{14}\\ {\mathit{s}}_{52}& {\mathit{s}}_{53}& {\mathit{s}}_{54}\\ {\mathit{s}}_{62}& {\mathit{s}}_{63}& {\mathit{s}}_{64}\end{array}\right]\left\{\begin{array}{c}{\overline{u}}_0\\ {\overline{v}}_0\\ {\overline{w}}_0\end{array}\right\}=\left\{\begin{array}{c}{\mathit{C}}_{\mathit{a}}\left(1\right)\\ {\mathit{C}}_{\mathit{a}}\left(5\right)\\ {\mathit{C}}_{\mathit{a}}\left(6\right)\end{array}\right\}.$$(33) , displacements, $\left[\begin{array}{ccc}{\overline{u}}_0& {\overline{v}}_0& {\overline{w}}_0\end{array}\right],$ at the lower surface in the sensor can be obtained; then using Eqs. (20)(20) $$\overline{\delta }\left(\overline{z}\right){=\mathit{e}}^{{\overline{G}}_{\mathit{k}}\overline{z}}{\mathit{e}}^{{-\overline{G}}_{\mathit{k}}{\overline{z}}_{\mathit{k}-1}}{\overline{\delta }}_{0\mathit{k}},\hspace{1em}{\overline{z}}_{\mathit{k}-1}\le \overline{z}\le {\overline{z}}_{\mathit{k}},$$(20) , (26-a) and (26-b) leads to determine mechanical and electrical state variables in three dimensions. Finally, to invert Laplace transforms of the states, modified Dubner’s and Abats’s numerical techniques are utilized to represent displacement and stress fields in the time domain (Abate and Whitt 1995).

Numerical results

A simply-supported square FGM plate with viscoelastic behavior sandwiched between two piezoelectric layers is considered. The material properties of the two piezoelectric layers are given in Table 1 (Alibeigloo 2010). The relaxation time constant is $\mathit{\tau }=5\mathit{s},$ the Poison ratio is $\mathit{v}=0.3$ and the parameters for the Prony series, related to Eq. (7)(7) $$\mathit{E}\left(\mathit{z},\mathit{t}\right){=\mathit{E}}_{1\mathit{m}}{\left(\frac{\mathit{z}}{{\mathit{h}}_{\mathit{p}}}\right)}^{{\mathit{m}}_1}{+\mathit{E}}_{2\mathit{m}}{\left(\frac{\mathit{z}}{{\mathit{h}}_{\mathit{p}}}\right)}^{{\mathit{m}}_2}{\mathit{e}}^{\left(-\frac{\mathit{t}}{\tau }\right)},$$(7) , as in Norouzi and Alibeigloo (2018) are $$E_{1m}=70\mathit{GPa},\hspace{1em}{E}_{2m}=105\mathit{GPa},\hspace{1em}{E}_{1c}=380\mathit{GPa},\hspace{1em}{E}_{2c}=570\mathit{GPa}$$

Table 1. Material properties of piezoelectric sensor and actuator (Alibeigloo 2010).

Elasticity moduli [10^9N m^−2]	$Q_{11}$	$Q_{12}$	$Q_{13}$	$Q_{22}$	$Q_{23}$	$Q_{33}$	$Q_{44}$	$Q_{55}$	$Q_{66}$
(PZT-4)	139	78	74	139	74	115	25.6	25.6	30.5
$(B{a}_{2}NaN{b}_5{O}_{15})$	239	104	50	274	52	135	65	66	76
Piezoelectric coefficients [cm^−2]	$e_{31}$	$e_{32}$	$e_{33}$	$e_{24}$	$e_{15}$
(PZT-4)	−5.2	−5.2	15.1	12.7	12.7
$(B{a}_{2}NaN{b}_5{O}_{15})$	−0.4	−0.3	4.3	3.4	2.8
Permittivities [10^−9 As V^−1m^−1]	${\epsilon }_{11}$	${\epsilon }_{22}$	${\epsilon }_{33}$
(PZT-4)	6.5	6.5	5.6
$(B{a}_{2}NaN{b}_5{O}_{15})$	1.96	2.01	0.28

It is noted that all presented results are their maximum value which are computed at $\mathit{t}=2\mathit{s}.$ The values of the other parameters involved in the plate are as follows unless the effects of variation of those parameters are studies. $$\mathit{a}/\mathit{b}=1,\hspace{1em}\mathit{a}/\mathit{h}=10,\hspace{1em}{\mathit{h}/\mathit{h}}_{\mathit{p}}=25,\mathit{\tau }=5\mathit{s}.$$

First, we study the convergence of the solution. We consider the length-to-thickness ratio as $a/h=10$ and the total thickness to piezoelectric layer thickness ratio as $h/h_p=25.$ The actuator is operated in a short circuit condition and only the uniform pressure ${\mathit{P}}_0$ is applied. The through-the-thickness distribution of the non-dimensional deflection, transverse normal stress, and transverse shear stress in the center of the plate are shown in Fig. 2(a)–(c) for various half wave numbers ($\mathit{m},\mathit{n}$) with 15 fictitious layers for the FGM core. From the figures, one can see that increasing the half wave numbers up to $\mathit{n}=\mathit{m}=25$ increases the transverse normal stress and the deflection while it decreases the transverse shear stress; all three entities are converged, which indicates fast convergence with respect to the half wave numbers ($\mathit{m},\mathit{n}$). Next, we keep $\mathit{n}=\mathit{m}=25$ fixed and study the convergence with respect to the number of fictitious layers for the FGM core.

Figure 2. Distribution of transverse normal stress, shear stress, and deflection along the thickness direction for different half wave numbers (m, n), using 15 fictitious layers for the FGM core layer.

Next, for validation of the present approach, an elastic square FGM plate under mechanical load, which is utilized in Brischetto et al. (2008), is adopted. Dimensionless stress and displacement components in the in-plane center of the plate at the bottom ($\mathit{z}=0$), top ($\mathit{z}=\mathit{h}$) and mid surface ($\mathit{z}=\mathit{h}/2$) of the plate is computed and compared with those reported in Brischetto et al. (2008). We use $\mathit{n}=\mathit{m}=25$ and $10$ fictitious layers. The results are presented in Table 2. According to Table 2, the present results are in good agreement with the results reported in Brischetto et al. (2008). The minor discrepancy is due to the CFU formulation which has been used in this reference.

Table 2. Comparison of stresses and displacements of an elastic square FGM plate under a mechanical load.

	$\frac{a}{h}=4$		$\frac{a}{h}=10$
	(Brischetto et al. 2008)	Analytical	(Brischetto et al. 2008)	Analytical
${\overline{w}}_{(h)}$	−13.46	−13.38	−168.9	−168.8
${\overline{w}}_{(\frac{h}{2})}$	−13.70	−13.77	−170.7	−170.39
${\overline{w}}_{(0)}$	−12.73	−12.80	−168.5	−168.3
${\overline{u}}_{(h)}$	4.02	4.09	26.17	26.165
${\overline{u}}_{(\frac{h}{2})}$	−0.08998	−0.090	0.7108	0.7121
${\overline{u}}_{(0)}$	−4.069	−4.08	−24.72	−24.78

The convergence of through-thickness distribution of transverse normal and shear stresses and deflection by increase in the number of fictitious layers is computed and depicted in Fig. 3(a)–(c). From the figures, it can be observed that the effect of increase in the number of mathematical layers on deflection is more significant. Moreover, it is seen that with 15 fictitious layers, all entities converge to constant values. It is noted that all of the upcoming results are based on this value for the mathematical layers.

Figure 3. Through-the-thickness distribution of the non-dimensional transverse normal stress, inplane normal stress, and transverse shear stress for different numbers of mathematical layers for the FGM core.

In the following results, we use $\mathit{n}=\mathit{m}=25$ and $15$ fictitious layers. The effect of the length-to-thickness ratio $\mathit{a}/\mathit{h}$ on the through-the-thickness distribution of dimensionless stresses and longitudinal displacement of the viscoelastic FGM plate are depicted in Fig. 4(a)–(c). One can see that an increase in the$\mathit{a}/\mathit{h}$ value causes the stiffness of the plate to decrease and consequently the stresses and the displacement to increase. Moreover, it is observed that continuity of the in-plane displacement and transverse shear stress at the interface of the FGM core and piezoelectric layer is satisfied. As expected, continuity of the in-plane shear stress at the interfaces is not satisfied. Additionally, it can be seen that the effect of variation in $\mathit{a}/\mathit{h}$ on ${\overline{\tau }}_{\mathit{\text{xy}}}$ is more significant near the top surface of the plate, whereas for the in-plane displacement this is the case near the bottom surface of the plate. Since the neutral surface in FGM core does not coincide with the mid surface, the maximum value of transverse shear stress is not located in the thickness center of the plate but shifted toward the top.

Figure 4. Effect of length to thickness ratio, a/h, on through-thickness distribution of nondimensional in-plane and transverse shear stresses and longitudinal displacement.

Effect of piezoelectric layer thickness on stress and displacement components in the absence of applied electric voltage are depicted in Fig. 5(a)–(c). According to the figures, in the absence of electric voltage, decreasing the piezoelectric layers thickness ratio up to ${\mathit{h}/\mathit{h}}_{\mathit{p}}=150$ causes the mechanical behavior of the viscoelastic FGM plate to converge to a constant value.

Figure 5. Effect of piezoelectric thick on distribution of non-dimensional in-plane and transverse shear stresses and longitudinal displacement.

Figure 6(a)–(d) represents the distribution of stress, displacements, and induced electric voltage of viscoelastic FGM plate with different relaxation time constant,$\mathit{\tau }.$ According to the figures, the displacements and electric voltage in the elastic plate are greater than those for the viscoelastic plate. Moreover, it is seen that in elastic FGM plate transverse shear stress near the mid surface is smaller and near the bottom and top surfaces is greater than that for the viscoelastic plate. Also, it can be observed that the increase causes the displacements and electric voltage to decrease.

Figure 6. Effect of time constant variation, τ, on the through-thickness distribution of nondimensional transverse shear stress, longitudinal and transverse displacements, and induced voltage.

Effect of applied voltage,${\overline{V}}_0,$ on the through-thickness distribution of longitudinal and transverse displacements of the square plate is depicted in Fig. 7(a)–(d). According to Fig. 7(a), voltage increase causes longitudinal displacement in both FGM and piezoelectric layers to increase linearly. From Fig. 7(b), it is seen that the increased voltage causes the deflection nonlinearly in the FGM layer to increase whereas it is linear in piezoelectric. Additionally, it is concluded that effect of applied voltage on deflection is more significant than that for the longitudinal displacement.

Figure 7. Influence applied voltage variation on the through-thickness distribution of nondimensional longitudinal and transverse displacements, in-plane normal stress, and transverse shear stress.

Displacement and stress distributions versus time with different relaxation time constants are shown in Fig. 8(a)–(c). According to the figures and as expected, stresses and displacement convergence toward a specified value. Additionally, it is concluded that an increase in the relaxation time constant causes the convergence time to increase. Also, it is seen that increasing relaxation time constant at each time causes the displacements to decrease and the stress to increase.

Figure 8. Time history of longitudinal displacement, transverse normal and shear stresses at the height of $\frac{(h_p+h_f+h_p)}{2}$ for different time constants.

Table 3 depicts dimensionless stresses and displacements for a thick, moderately thick and thin viscoelastic rectangular FGM plate embedded with piezoelectric layers, with different relaxation time constants. Regardless of relaxation time constant, increase in the length to thickness ratio, $\mathit{a}/\mathit{h},$ causes an increase in the displacement and stress components. From Table 3, it can be seen that increase relaxation time constant causes to increase stresses whereas displacements decrease. Moreover, it is concluded that effect of relaxation time constant in thin plate is more significant.

Table 3. Maximum stresses and displacements of the sandwich FGM plate in elastic ($\mathrm{\tau }=0$) and viscoelastic ($\mathrm{\tau }\mathrm{>}0$) cases for different $\frac{a}{h}$ ratio.

$\frac{a}{h}$	$\tau$	${\overline{\sigma }}_x$	$\overline{u}$	$\overline{w}$	${\overline{\tau }}_{xz}$	${\overline{\tau }}_{xy}$
10	0 (s)	33.53	0.3	1.7	5.0	−22.55
	5 (s)	36.62	0.2	0.8	5.1	−24.68
	10 (s)	36.76	0.2	0.8	5.1	−24.77
20	0 (s)	133.78	2.7	25.7	10	−90.39
	5 (s)	143.13	1.3	12.5	10.3	−98.78
	10 (s)	146.70	1.2	11.9	10.3	−99.17
40	0 (s)	534.91	21.3	407.1	20.0	−362.05
	5 (s)	584.31	10.4	198	20.6	−395.54
	10 (s)	586.57	9.9	188	20.6	−397.10

Conclusion

In the present paper, a three-dimensional static analysis of an FGM viscoelastic plate sandwiched between two piezoelectric layers and subjected to electromechanical loading was conducted within the framework of linearized three-dimensional viscoelasticity. The analytical solution was computed by applying a state-space approach through the thickness direction and a Fourier series approach for the in-plane problem. The governing state-space equations were solved in the Laplace domain and the response in the time domain was obtained using the inverse Laplace transform. Based on the numerical results, it was demonstrated that the viscoelastic properties strongly affect the bending behavior of the FGM plate and the following conclusions were derived:

• Mechanical displacements and electric voltage in the purely elastic plate are larger whereas the stresses are smaller than those for the viscoelastic plate.

• Increase of the relaxation time constant causes the displacements, the induced electric voltage, and the stresses to increase.

• The rate of convergence to elastic behavior decreases with increasing relaxation time constant.

• In the viscoelastic plate, the effect of the relaxation time constant increases with increasing length-to-thickness ratio.

• The effect of the applied voltage on the deflection compared to the in-plane displacements is more significant.

• The through-the-thickness distribution of the deflection due to an applied voltage is nonlinear whereas it is linear in the piezoelectric layers.

• The effect of increasing length-to-thickness ratio on the behavior of the viscoelastic plate more pronounced than that for the elastic plate.

Nomenclature
a, b, h	=	plate dimensions in x-, y- and z-directions
${\mathit{C}}_{\mathit{\text{ij}}}(\mathit{i},\mathit{j}=1.2,\dots .,6)$	=	relaxation moduli coefficients
${\mathit{D}}_{\mathit{x}}{,\mathit{D}}_{\mathit{y}}{,\mathit{D}}_{\mathit{z}}$	=	electric displacement
${\mathit{E}}_{\mathit{x}}{,\mathit{E}}_{\mathit{y}}{,\mathit{E}}_{\mathit{z}}$	=	electric field in x-, y- and z-directions
$\mathit{E}$	=	Young’s modulus
$\mathit{e},\mathit{\eta }$	=	piezoelectric constants and dielectric constants
${\mathit{d}}_1$	=	piezoelectric modulus
${\mathit{h}}_{\mathit{f}}{,\mathit{h}}_{\mathit{p}}$	=	thicknesses of the FGM and piezoelectric layers
$\mathit{n},\mathit{m}$	=	half wave numbers in the x-, y- and z-directions
$\mathit{u},\mathit{v},$ and $\mathit{w}$	=	displacement components in the x-, y- and z-directions
${\sigma }_x,{\sigma }_y,{\sigma }_z$	=	normal stresses
${\tau }_{xy},{\tau }_{zx},{\tau }_{zy}$	=	shear stresses
${\epsilon }_x,{\epsilon }_y,{\epsilon }_z$	=	normal strains
${\gamma }_{xy},{\gamma }_{zx},{\gamma }_{zy}$	=	shear strains
$\mathit{\tau }$	=	relaxation time constant
${\mathit{\delta }}_{\mathit{f}},$ ${\mathit{\delta }}_{\mathit{p}}$	=	state vectors of the FGM and piezoelectric layers
$\mathit{\psi }$	=	electric voltage
$\mathit{v}$	=	Poison ratio

Acknowledgment

The authors acknowledge the TU Wien Bibliothek for financial support through its Open Access Funding Program.

Appendix

In the present paper the piezoelectric layers are assumed to be orthorhombic piezoelectric materials of crystal class 2 mm. (A-1) $$\mathbf{C}=\left[\begin{array}{cccccc}\frac{\left(1-\nu \right)E\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& 0& 0& 0\\ \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{\left(1-\nu \right)E\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& 0& 0& 0\\ \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{\left(1-\nu \right)E\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& 0& 0& 0\\ 0& 0& 0& \frac{E\left(z,t\right)}{2\left(1+\nu \right)\left(1-2\nu \right)}& 0& 0\\ 0& 0& 0& 0& \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}& 0\\ 0& 0& 0& 0& 0& \frac{vE\left(z,t\right)}{\left(1+\nu \right)\left(1-2\nu \right)}\end{array}\right]$$(A-1) (A-2) $$G_f=\left[\begin{array}{cccccc}0& 0& 0& 0& -\frac{\partial }{\partial x}& -\frac{\partial }{\partial y}\\ 0& 0& 0& -\frac{\partial }{\partial x}& \frac{2\left(1+\nu \right)}{s\widehat{E}}& 0\\ 0& 0& 0& -\frac{\partial }{\partial y}& 0& \frac{2\left(1+\nu \right)}{s\widehat{E}}\\ \frac{\left(1+\nu \right)\left(1-2\nu \right)}{\left(1-\nu \right)s\widehat{E}}& -\frac{\nu }{\left(1-\nu \right)}\frac{\partial }{\partial x}& -\frac{\nu }{\left(1-\nu \right)}\frac{\partial }{\partial y}& 0& 0& 0\\ \frac{-\nu }{\left(1-\nu \right)}\frac{\partial }{\partial x}& \frac{-s\widehat{E}}{1-{\nu }^2}\frac{{\partial }^2}{\partial x^2}+\frac{s\widehat{E}}{2\left(1+\nu \right)}\frac{{\partial }^2}{\partial y^2}& -\frac{s\widehat{E}}{2\left(1-\nu \right)}\frac{\partial }{\partial y\partial x}& 0& 0& 0\\ \frac{-\nu }{\left(1-\nu \right)}\frac{\partial }{\partial y}& -\frac{s\widehat{E}}{2\left(1-\nu \right)}\frac{\partial }{\partial y\partial x}& \frac{-s\widehat{E}}{1-{\nu }^2}\frac{{\partial }^2}{\partial y^2}+\frac{s\widehat{E}}{2\left(1+\nu \right)}\frac{{\partial }^2}{\partial x^2}& 0& 0& 0\end{array}\right]$$(A-2) (A-3) $$C_p=\left[\begin{array}{ccc}\begin{array}{cc}{Q}_{11}& Q_{12}\\ Q_{12}& Q_{22}\end{array}& \begin{array}{cc}{Q}_{13}& 0\\ Q_{23}& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}{Q}_{13}& Q_{23}\\ 0& 0\end{array}& \begin{array}{cc}{Q}_{33}& 0\\ 0& Q_{44}\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}{Q}_{55}& 0\\ 0& Q_{66}\end{array}\end{array}\right]\mathit{}\mathit{},\mathit{}\mathit{}{\epsilon }_p=\left\{\begin{array}{c}{\widehat{\epsilon }}_{\mathit{x}}\\ {\widehat{\epsilon }}_y\\ {\widehat{\epsilon }}_z\\ {\widehat{\gamma }}_{xz}\\ {\widehat{\gamma }}_{yz}\\ {\widehat{\gamma }}_{xy}\end{array}\right\}\mathit{},\mathit{}\mathit{}{\sigma }_p=\left\{\begin{array}{c}{\widehat{\sigma }}_x\\ {\widehat{\sigma }}_y\\ {\widehat{\sigma }}_z\\ {\widehat{\tau }}_{xz}\\ {\widehat{\tau }}_{yz}\\ {\widehat{\tau }}_{xy}\end{array}\right\}$$(A-3) (A-4) $$D={\left\{\begin{array}{ccc}{D}_x& D_y& D_z\end{array}\right\}}^T{\mathit{},\mathit{\text{E}}}_p={\left\{\begin{array}{ccc}{E}_x& E_y& E_z\end{array}\right\}}^T,\mathit{}\eta =\left[\begin{array}{ccc}{\eta }_1& 0& 0\\ 0& {\eta }_2& 0\\ 0& 0& {\eta }_3\end{array}\right]$$(A-4) (A-5) $$e=\left[\begin{array}{c}\begin{array}{cccc}0& 0& 0& \begin{array}{ccc}0& e_5& 0\end{array}\end{array}\\ \begin{array}{cccc}0& 0& 0& \begin{array}{ccc}{e}_4& 0& 0\end{array}\end{array}\\ \begin{array}{cccc}{e}_1& e_2& e_3& \begin{array}{ccc}0& 0& 0\end{array}\end{array}\end{array}\right]$$(A-5) (A-6) $$G_P=\left[\begin{array}{cccccccc}0& 0& 0& 0& -\frac{\partial }{\partial x}& -\frac{\partial }{\partial y}& 0& 0\\ 0& 0& 0& -\frac{\partial }{\partial x}& \frac{1}{Q_{55}}& 0& 0& -\frac{e_5}{Q_{55}}\frac{\partial }{\partial x}\\ 0& 0& 0& -\frac{\partial }{\partial y}& 0& \frac{1}{Q_{44}}& 0& -\frac{e_4}{Q_{44}}\frac{\partial }{\partial y}\\ \frac{{\eta }_3}{Q_{33}{\eta }_3+e_3^2}& -\frac{Q_{13}{\eta }_3+e_1{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial x}& G_{P43}\frac{\partial }{\partial y}& 0& 0& 0& G_{P47}& 0\\ -\frac{Q_{13}{\eta }_3+e_1{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial x}& G_{P52}& G_{P53}& 0& 0& 0& \frac{Q_{33}{e}_1-C_{13}{e}_3}{Q_{33}{\eta }_3+e_3^2}& 0\\ -\frac{Q_{23}{\eta }_3+e_2{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial y}& G_{P62}& G_{P63}& 0& 0& 0& \frac{Q_{33}{e}_1-Q_{23}{e}_3}{Q_{33}{\eta }_3+e_3^2}\frac{\partial }{\partial y}& 0\\ 0& 0& 0& 0& -\frac{e_5}{Q_{55}}\frac{\partial }{\partial x}& -\frac{e_4}{Q_{44}}\frac{\partial }{\partial y}& 0& \left(\frac{e_5^2}{Q_{55}}+{\eta }_1\right)\frac{{\partial }^2}{\partial x^2}+\left(\frac{e_4^2}{Q_{44}}+{\eta }_2\right)\frac{{\partial }^2}{\partial y^2}\\ \frac{e_3}{Q_{33}{\eta }_3+e_3^2}& -\left(\frac{Q_{13}{e}_3-Q_{33}{e}_1}{Q_{33}{\eta }_3+e_3^2}\right)\frac{\partial }{\partial x}& -\left(\frac{Q_{23}{e}_3-Q_{33}{e}_2}{Q_{33}{\eta }_3+e_3^2}\right)\frac{\partial }{\partial y}& 0& 0& 0& -\frac{Q_{33}}{Q_{33}{\eta }_3+e_3^2}& 0\end{array}\right]$$(A-6) $$\begin{array}{c}{G}_{P43}=-\frac{Q_{23}{\eta }_3+e_2{e}_3}{Q_{33}{\eta }_3+e_3^2},{\mathit{\text{G}}}_{P47}=\frac{e_3}{Q_{33}{\eta }_3+e_3^2},{\mathit{\text{G}}}_{P52}=\left(-Q_{11}+\frac{Q_{13}^2{\eta }_3+e_1{e}_3{Q}_{13}+e_1{e}_3{Q}_{13}-Q_{33}{e}_1^2}{Q_{33}{\eta }_3+e_3^2}\right)\frac{{\partial }^2}{\partial x^2}-Q_{66}\frac{{\partial }^2}{\partial y^2}\\ G_{P62}=G_{P53}=\left(-Q_{12}-Q_{66}+\frac{Q_{13}{Q}_{23}{\eta }_3+e_2{e}_3{Q}_{13}+e_1{e}_3{Q}_{23}-e_1{e}_2{Q}_{33}}{Q_{33}{\eta }_3+e_3^2}\right)\frac{{\partial }^2}{\partial x\partial y}\\ G_{P63}=\left(-Q_{22}+\frac{Q_{23}^2{\eta }_3+e_2{e}_3{Q}_{23}+e_2{e}_3{Q}_{23}-Q_{33}{e}_2^2}{Q_{33}{\eta }_3+e_3^2}\right)\frac{{\partial }^2}{\partial y^2}\end{array}$$ (A-7) $${\overline{G}}_f=\left[\begin{array}{cccccc}0& 0& 0& 0& \left(\frac{h}{a}\right){\overline{P}}_m& \left(\frac{h}{b}\right){\overline{P}}_n\\ 0& 0& 0& -\left(\frac{h}{a}\right){\overline{P}}_m& \frac{2\left(1+\nu \right)}{s\overline{E}}& 0\\ 0& 0& 0& -\left(\frac{h}{b}\right){\overline{P}}_n& 0& \frac{2\left(1+\nu \right)}{s\overline{E}}\\ \frac{\left(1+\nu \right)\left(1-2\nu \right)}{\left(1-\nu \right)s\overline{E}}& {\overline{P}}_m\left(\frac{h}{a}\right)\frac{\nu }{\left(1-\nu \right)}& {\overline{P}}_n\left(\frac{h}{b}\right)\frac{\nu }{\left(1-\nu \right)}& 0& 0& 0\\ -{\overline{P}}_m\left(\frac{h}{a}\right)\frac{\nu }{\left(1-\nu \right)}& {\overline{P}}^2{}_m{\left(\frac{h}{a}\right)}^2\frac{s\overline{E}}{\left(1-{\nu }^2\right)}-{\overline{P}}^2{}_n{\left(\frac{h}{b}\right)}^2\frac{s\overline{E}}{2\left(1+\nu \right)}& {\overline{P}}_m{\overline{P}}_n\left(\frac{h}{a}\right)\left(\frac{h}{b}\right)\frac{s\overline{E}}{2\left(1+\nu \right)}& 0& 0& 0\\ -{\overline{P}}_n\left(\frac{h}{b}\right)\frac{\nu }{\left(1-\nu \right)}& {\overline{P}}_m{\overline{P}}_n\left(\frac{h}{a}\right)\left(\frac{h}{b}\right)\frac{s\overline{E}}{2\left(1+\nu \right)}& {\overline{P}}^2{}_n{\left(\frac{h}{b}\right)}^2\frac{s\overline{E}}{\left(1-{\nu }^2\right)}+{\overline{P}}^2{}_m{\left(\frac{h}{a}\right)}^2\frac{s\overline{E}}{2\left(1+\nu \right)}& 0& 0& 0\end{array}\right]$$(A-7) (A-8) $$\overline{F}=\left[\begin{array}{ccc}\frac{\nu }{1-\nu }& \left(\frac{h}{a}\right){\overline{P}}_m\frac{s\overline{E}}{1-{\nu }^2}& \left(\frac{h}{b}\right){\overline{P}}_n\frac{s\overline{E}}{1-{\nu }^2}\\ \frac{\nu }{1-\nu }& \left(\frac{h}{a}\right){\overline{P}}_m\frac{s\overline{E}}{1-{\nu }^2}& \left(\frac{h}{b}\right){\overline{P}}_n\frac{s\overline{E}}{1-{\nu }^2}\\ 0& \left(\frac{h}{b}\right){\overline{P}}_n\frac{s\overline{E}}{2(1+\nu )}& \left(\frac{h}{a}\right){\overline{P}}_m\frac{s\overline{E}}{2(1+\nu )}\end{array}\right]$$(A-8) $${\overline{G}}_P=\left[\begin{array}{cccccccc}0& 0& 0& 0& \left(\frac{h}{a}\right){\overline{P}}_m& \left(\frac{h}{b}\right){\overline{P}}_n& 0& 0\\ 0& 0& 0& -\left(\frac{h}{a}\right){\overline{P}}_m& \frac{1}{{\overline{Q}}_{55}}& 0& 0& -\frac{{\overline{e}}_5}{{\overline{Q}}_{55}}{\overline{P}}_m\\ 0& 0& 0& -\left(\frac{h}{b}\right){\overline{P}}_n& 0& \frac{1}{{\overline{Q}}_{44}}& 0& -\frac{{\overline{e}}_4}{{\overline{Q}}_{44}}{\overline{P}}_n\\ \frac{{\overline{\eta }}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}& \frac{{\overline{Q}}_{13}{\overline{\eta }}_3+{\overline{e}}_1{\overline{e}}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\left(\frac{h}{a}\right){\overline{P}}_m& {\overline{G}}_{P43}\left(\frac{h}{b}\right){\overline{P}}_n& 0& 0& 0& {\overline{G}}_{P47}& 0\\ -\frac{{\overline{Q}}_{13}{\overline{\eta }}_3+{\overline{e}}_1{\overline{e}}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\left(\frac{h}{a}\right){\overline{P}}_m& {\overline{G}}_{P52}& {\overline{G}}_{P53}& 0& 0& 0& \frac{{\overline{Q}}_{33}{\overline{e}}_1}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}{\overline{P}}_m& 0\\ -\frac{{\overline{Q}}_{23}{\overline{\eta }}_3+{\overline{e}}_2{\overline{e}}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}{\overline{P}}_n& {\overline{G}}_{P62}& {\overline{G}}_{P63}& 0& 0& 0& \frac{{\overline{e}}_2{\overline{Q}}_{33}-{\overline{e}}_3{\overline{Q}}_{23}}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}{\overline{P}}_n& 0\\ 0& 0& 0& 0& \frac{{\overline{e}}_5}{{\overline{Q}}_{55}}\left(\frac{h}{a}\right){\overline{P}}_m& \frac{{\overline{e}}_4}{{\overline{Q}}_{44}}\left(\frac{h}{b}\right){\overline{P}}_n& 0& -\left(\frac{{\overline{e}}_5^2}{{\overline{Q}}_{55}}+{\overline{\eta }}_1\right){\overline{P}}^2{}_m-\left(\frac{{\overline{e}}_4^2}{{\overline{Q}}_{44}}+{\overline{\eta }}_2\right){\overline{P}}^2{}_n\\ \frac{{\overline{e}}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}& \left(\frac{{\overline{Q}}_{13}{\overline{e}}_3-{\overline{Q}}_{33}{\overline{e}}_1}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right)\left(\frac{h}{a}\right){\overline{P}}_m& \left(\frac{{\overline{Q}}_{23}{\overline{e}}_3-{\overline{Q}}_{33}{\overline{e}}_2}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right)\left(\frac{h}{b}\right){\overline{P}}_n& 0& 0& 0& -\frac{{\overline{Q}}_{33}}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}& 0\end{array}\right]$$ (A-9) $$\begin{array}{c}{\overline{G}}_{P43}=\frac{{\overline{Q}}_{23}{\overline{\eta }}_3+{\overline{e}}_2{\overline{e}}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2},{\overline{G}}_{P47}=\frac{e_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2},\\ {\overline{G}}_{P52}=\left({\overline{Q}}_{11}-\frac{{\overline{Q}}_{13}^2{\overline{\eta }}_3+{\overline{e}}_1{\overline{e}}_3{\overline{Q}}_{13}+{\overline{e}}_1{\overline{e}}_3{\overline{Q}}_{13}-{\overline{Q}}_{33}{\overline{e}}_1^2}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right){\left(\frac{h}{a}\right)}^2{\overline{P}}^2{}_m+{\overline{Q}}_{66}{\left(\frac{h}{b}\right)}^2{\overline{P}}^2{}_n\\ {\overline{G}}_{P62}={\overline{G}}_{P53}=\left(\frac{h}{a}\right)\left(\frac{h}{b}\right)\left({\overline{Q}}_{12}+{\overline{Q}}_{66}-\frac{{\overline{Q}}_{13}{\overline{Q}}_{23}{\overline{\eta }}_3+{\overline{e}}_2{\overline{e}}_3{\overline{Q}}_{13}+{\overline{e}}_1{\overline{e}}_3{\overline{Q}}_{23}-{\overline{e}}_1{\overline{e}}_2{\overline{Q}}_{33}}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right){\overline{P}}_m{\overline{P}}_n\\ {\overline{G}}_{P63}={\overline{Q}}_{66}{\left(\frac{h}{a}\right)}^2{\overline{P}}_m^2+\left({\overline{Q}}_{22}-\frac{{\overline{Q}}_{23}^2{\overline{\eta }}_3+{\overline{e}}_2{\overline{e}}_3{\overline{Q}}_{23}+{\overline{e}}_2{\overline{e}}_3{\overline{Q}}_{23}-{\overline{Q}}_{33}{\overline{e}}_2^2}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right){\left(\frac{h}{b}\right)}^2{\overline{P}}_n^2\end{array}$$(A-9) (A-10) $$\begin{array}{c}{\overline{F}}_p={\left[F_{ij}\right]}_{4\times 4},{\mathit{\text{F}}}_{11}=\frac{{\overline{Q}}_{13}{\overline{\eta }}_3+{\overline{e}}_1{\overline{e}}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2},\\ F_{12}=\left(\frac{h}{a}\right){\overline{P}}_m\left(-{\overline{Q}}_{11}+\frac{{\overline{Q}}_{13}({\overline{Q}}_{13}{\overline{\eta }}_3+{\overline{e}}_1{\overline{e}}_3)+{\overline{e}}_1({\overline{e}}_3{\overline{Q}}_{13}-{\overline{e}}_1{\overline{Q}}_{33})}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right),\\ F_{13}=\left(\frac{h}{b}\right){\overline{P}}_n\left(-{\overline{Q}}_{12}+\frac{{\overline{Q}}_{13}({\overline{Q}}_{23}{\overline{\eta }}_3+{\overline{e}}_2{\overline{e}}_3)+{\overline{e}}_1({\overline{e}}_3{\overline{Q}}_{23}-{\overline{e}}_2{\overline{Q}}_{33})}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right),\\ F_{14}=\frac{{\overline{Q}}_{13}{\overline{e}}_3+{\overline{e}}_1{\overline{Q}}_{33}}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2},{\mathit{\text{F}}}_{21}=\frac{{\overline{Q}}_{23}{\overline{\eta }}_3+{\overline{e}}_2{\overline{e}}_3}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2},\\ F_{22}=\left(\frac{h}{a}\right){\overline{P}}_m\left(-{\overline{Q}}_{12}+\frac{{\overline{Q}}_{23}({\overline{Q}}_{13}{\overline{\eta }}_3+{\overline{e}}_1{\overline{e}}_3)+{\overline{e}}_2({\overline{e}}_3{\overline{Q}}_{13}-{\overline{e}}_1{\overline{Q}}_{33})}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right),\\ F_{23}=\left(\frac{h}{b}\right){\overline{P}}_n\left(-{\overline{Q}}_{22}+\frac{{\overline{Q}}_{23}({\overline{Q}}_{23}{\overline{\eta }}_3+{\overline{e}}_2{\overline{e}}_3)+{\overline{e}}_2({\overline{e}}_3{\overline{Q}}_{23}-{\overline{e}}_2{\overline{Q}}_{33})}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2}\right),\\ F_{24}=\frac{{\overline{Q}}_{23}{\overline{e}}_3+{\overline{e}}_2{\overline{Q}}_{33}}{{\overline{Q}}_{33}{\overline{\eta }}_3+{\overline{e}}_3^2},{\mathit{\text{F}}}_{31}=0,F_{32}=\left(\frac{h}{a}\right){\overline{P}}_n{\overline{Q}}_{66},{\mathit{\text{F}}}_{33}=\left(\frac{h}{a}\right){\overline{P}}_m{\overline{Q}}_{66},k_{34}=0.\end{array}$$(A-10)