An Efficient Spectral Element Model with Electric DOFs for the Static and Dynamic Analysis of a Piezoelectric Bimorph

An efficient spectral element (SE) model for static and dynamic analysis of a piezoelectric bimorph is proposed. It combines an equivalent single layer (ESL) model for the mechanical displacement field with a sublayer approximation for the electric potential. The 2D Gauss-Lobatto-Legendre (GLL) shape functions are used to discretize the displacements and then the governing equation of motion is derived following the standard SE method procedure. It is shown numerically that the present SE model can well predict both the global and local responses such as mechanical displacements, natural frequencies, and the electric potentials across the bimorph thickness. In the case of bimorph sensor application, it is revealed that the distribution of the induced electric potential across the thickness does not affect the global natural frequencies much. Furthermore, the effects of the order of Legendre polynomial and the mesh size on the convergence rate are investigated. Comparison of the present results for a bimorph sensor with those from 3D finite element (FE) simulations establishes that the present SE model is accurate, robust, and computationally efficient.


Introduction
Piezoelectric materials, especially lead zirconate titanate (PZT), can function either as actuator or as sensor for their inherent coupling electromechanical character.Due to its wide ranging applications in electroacoustic transducers, medical devices, microrobot, and atomic force microscope (AFM) cantilevers, the study of smart structures consisting of PZT sensors and PZT actuators has drawn considerable attention of many researchers in the fields of active vibration control, noise attenuation, and damage detection.The most popular simple PZT sensor or actuator consists usually of a slab of piezoelectric ceramic such that the PZT layer expands or contracts mainly in its length direction.However, the motion of a single layer, as well as the induced electric potential, is extremely small.To achieve practically meaningful actuation or sensing capabilities in PZT devices, a piezoelectric bimorph consisting of two PZT layers is commonly used for the reason that it can produce flexural deformation significantly larger than the length or thickness deformation of the individual layers [1][2][3][4][5][6].However, the applications of the piezoelectric bimorph require the development of admissible approaches entailing capabilities to predict the global responses of the bimorph structure, such as the deflection and natural frequencies.Additionally, the approaches should address the local responses, such as the through-the-thickness variation of the electric potentials.
In the past three decades, a large variety of models have been developed to predict the static and dynamic responses of piezoelectric bimorph structures under all kinds of electromechanical loads with emphasis on approximating the mechanical displacement and electric potential correctly.The classification of the various models is mainly based on the kinematic assumption for approximating the through-thethickness variation of the electromechanical state variables and representation method of the piezoelectric layers [1].The accurate responses of the piezoelectric bimorph structures can be obtained by solving the 3-dimensional (3D) coupled field equations with exact satisfaction of the mechanical and electric boundary conditions [7].However, the 3D analytical solutions do not provide the results for more general case of loading and complicated boundary conditions.Thus, there is a genuine requirement of improved numerical techniques such as finite element (FE) method.The 3D FE modeling of piezoelectric bimorph structures, which is now available in many commercial FE software programs such as ABAQUS and ANSYS, will result in systems with a large number of degrees of freedom (DOFs).To overcome the computational inefficiency associated with the 3D FE models, the equivalent single layer (ESL) model applied to both mechanical and electrical unknowns has been proposed.There are two main kinds of theories used for ESL model.One is the classical laminated plate theory (CLPT) [8,9], and the other is the shear deformation theory, which branches out into first-order shear deformation theory (FSDT) [10,11] and higher-order shear deformation theory (HSDT) [12,13].The ESL model is simple and capable of predicting the global responses of the bimorph, but it does not account for the nonlinear distribution of the electric potential across piezoelectric layers as observed from 3D solutions [4].This shortcoming inspired the researchers to develop layer-wise theory [4,[14][15][16] or the sublayer theory [2,7,[17][18][19] for approximation of the electric potential.In the latter case, the piezoelectric layer is divided into appropriate number of sublayers along the thickness direction and a linear through-the-thickness electric potential distribution for each sublayer is assumed.It is expected that the actual distribution of the electric potential across the bimorph thickness can be approached with more sublayers adopted.Therefore, a sensible strategy is to use an ESL model for the mechanical variables and a layerwise theory or a sublayer theory for the electric variables, respectively [20][21][22].
The FE model which considers a global ESL approximation for the mechanical field variables and a layer-wise theory or a sublayer theory for the electric potential will inevitably result in a large number of DOFs for practical dynamical problems.Recently, the spectral element (SE) method, which combines the geometric flexibility of the FE method with the accuracy of the global pseudospectral method, is widely used in many research areas related to seismology, fluid dynamics, and acoustics [23][24][25].More recently, the SE method was used to simulate wave propagations for the purpose of damage detection in structures [26].In fact, the SE method and FE method are closely related and built on the same ideas.The main difference between them is that SE method uses orthogonal polynomials, such as Legendre and Chebyshev polynomials, in the approximation functions; therefore, the mass matrix is diagonal, which is a very significant advantage over conventional FE method especially for dynamic analysis.Moreover, to have an accurate simulation with conventional FE method, a mesh with large number of nodes and elements is inevitably needed.The SE method, in which the polynomial order is increased and the mesh size is decreased, can be used to overcome this problem.However, it appears that the use of SE method for problems of static and dynamic analysis of piezoelectric bimorph has not been widely reported in the literatures so far.
The objective of the present work is to develop an efficient and accurate electromechanical coupled 2D SE model for static and dynamic analysis of a piezoelectric bimorph using combination of ESL for the displacements and a sublayer model for the electric potentials.The remainder of this paper is outlined as follows.The approximations for the displacements and the electric potentials are given in Section 2; the SE model with electric DOFs is then derived utilizing Legendre orthogonal polynomials in the interpolation function.In Section 3, the global responses, such as the deflection and natural frequencies, and the local responses, such as the electric potential across the thickness, are presented.The results are also compared to those provided by the 3D FE simulations.The convergence of the present SE model and the robustness of algorithm with respect to the order of the Legendre polynomial as well as mesh size are then investigated.At last the closing remarks and discussion of the results are given in Section 4.

Mathematical Formulation
2.1.Constitutive Relationships, Displacement, and Strain.A piezoelectric bimorph made of two PZT-4 piezoelectric layers (as shown in Figure 1) undergoing a surface density of normal force and electric potential applied to the top and bottom surfaces is considered in this work.The length  and width  of the bimorph are 25 mm and 12.5 mm, respectively.Both piezoelectric layers have the same thickness 0.5ℎ and are poled in the -direction.The reference - plane is chosen to be the middle plane of the bimorph, and the -axis is defined as the direction normal to the plane according to the right-hand rule.Generally, the linear constitutive equations of piezoelectric materials, including the converse and direct piezoelectric effects, can be written as where T represent stress vector and strain vector, respectively.
T is the electric displacement vector, c is the elastic coefficient matrix, g is the dielectric coefficient matrix, and e is the piezoelectric stress coefficient matrix.
The FSDT is based on the constant transverse shear strain assumption, which leads to the displacement field [27]  (, , , ) =  (, , ) +  (, , ) , where , V, and  are the displacements in the -, -, and -directions, respectively., V, and  are the in-plane and transverse displacements of a point (, ) on the middle plane, respectively. and - denote the rotations of a  transverse normal about the and -axes, respectively.We define where U is the displacement vector and U is a generalized displacement vector.Then (2) can be expressed in matrix form as where The relationship between the strains and the displacements can be written as follows: where L is the derivation operator defined as where   is the th order Legendre polynomial.In fact, the nodes are the 2D Gauss-Lobatto-Legendre (GLL) points.In this way the nodes of the element can be specified in the local coordinate system of the element, as shown in Figure 2. The 2D shape function built on the specified node (  ,   ) can be written as where ℎ  () is the 1D shape function of order  at the 1D GLL points   which can be defined as [26] ℎ The 2D shape functions are used for interpolating both the element coordinates and the element displacements.Consequently, coordinates  and  within each Ω e may be uniquely related to  and  upon the invertible mapping where   and   denote the coordinates of  and , respectively, of the element nodes (  ,   ).The generalized displacements , V, , , and  over a reference element Ω ref are discretized by the 2D shape functions as ⟨ (, ) , V (, ) ,  (, ) ,  (, ) ,  (, )⟩ where   , V  ,   ,   , and   are the nodal values of the generalized displacements.The discrete element nodal displacement vector is expressed as where q , is the displacement vector of the node (  ,   ) Substituting ( 12) into (4) yields where N  is the matrix of displacement shape function which can be written as with where I 5×5 is a 5 × 5 identity matrix.Substituting ( 15) into (6) yields where B  is strain-displacement matrix which can be written as

A Sublayer Model for the Electric Potentials.
In order to accurately model the through-the-thickness distribution of the electric potential, each layer of the piezoelectric bimorph is subdivided mathematically into  thinner sublayers.As shown in Figure 1, the sublayers are numbered in top-tobottom order.The -coordinates of the top and bottom surfaces of the th sublayer are denoted by   and  −1 , respectively.It is assumed that in each sublayer the electric potential   () has a linear variation across the thickness such that where N   is the shape function of the electric potential and Φ  is a column matrix composed of the electric potentials at the top and the bottom surfaces of the th sublayer, which can be expressed as In this way the electric potential is approximated as piecewise linear across the thickness and it is expected that the actual nonlinear electric potential field of the piezoelectric bimorph can be approached with an appropriate number of sublayers.
Considering that the piezoelectric bimorph is discretized using 2D mesh, each sublayer is also discretized using the same mesh to keep the compatibility.Consequently, an element potential vector Φ e which is of the following form is then introduced in the spectral plate finite element: The surface potential of the sublayer,   , is assumed to be constant over the element and  0 ,  1 , . . .,  2 are treated as elemental degrees of freedom (DOFs), as illustrated in Figure 2. Furthermore, the top and bottom surfaces of the piezoelectric layers are always coated with metallic coatings of zero thickness and the potentials on the electrodes should be taken as independent of , .Thus the present approach combines an ESL theory for the displacement field and a piecewise linear approximation for the electric potential.
Under the quasi electrostatic approximation, the electric field and the electric potential in each sublayer have the following relationship: where E  () is the electric field of the th sublayer and B   is the electric field-potential matrix, given by 2.4.Governing Equations.Hamilton's variational principle is adopted in the derivation of the elementary governing equation of motion q e + K e  q e + K e  Φ e = F e  , where M e  denotes the elementary mass matrix; K e  is mechanical stiffness matrix; K e  and K e  are the piezoelectric coupling matrices; K e  is the dielectric permittivity matrix; F e  and F e  denote the nodal external force vector and the nodal externally applied charge vector, respectively: where  is the mass density, P  is the surface force vector, and q  is the surface charge density vector.J is the Jacobian matrix of the mapping (11) which is defined by The GLL integration rule is then used to calculate the characteristic matrices and the nodal force vector in ( 25) at the elemental level [26].
In this study, the interface between the two PZT layers is grounded.The bimorph is supposed to be used as a sensor in closed circuit, in which the top and bottom surfaces are grounded and a uniform pressure load is applied to the upper surface.By applying the electric boundary conditions, the DOFs for the electric potentials are condensed out such that ( 25) is finally of the form where K e  is the mechanical stiffness matrix induced by the electromechanical coupling of the piezoelectric materials and F e  is the mechanical forces induced by the applied voltages of piezoelectric actuators [2].The electric potential is then recovered by the inverse process of the aforementioned condensation.Assembling all elementary equations, one can have a global dynamic system equation where M  , K  , K  , F  , and F  are the assembled counterparts of matrices M e  , K e  , K e  , F e  , and F e  ; q is the global nodal displacement vector.Since the electric potential DOFs for the sublayers have been condensed out, this approach will not result in an excessive number of potential field variables.For the purpose of static analysis, the governing equations of motion in (29) reduces to For the dynamic frequency analysis problems, the corresponding eigenvalue problem is where  is the natural frequency.

Numerical Results and Discussion
In this section, case studies are carried out to demonstrate the efficiency and accuracy of the present model in estimating both the global responses, such as the deflection and natural frequencies, and the local responses, such as the through-the-thickness variation of the electric potential.A simply supported rectangular piezoelectric bimorph shown in Figure 1, which was analyzed by Fernandes and Pouget [1], is considered here.The material constants of PZT-4 are given as Different slenderness ratios such as  = /ℎ = 5,  = 10, and  = 50, which represent the thick, moderately thick, and thin bimorph plate, respectively, are considered.Unless otherwise stated, the order of Legendre polynomial is chosen as  = 5, and the mesh in Figure 2 is used in this work.

Static Responses of the Bimorph
Sensor.The bimorph is supposed to be used as a sensor in closed circuit.To achieve practically meaningful sensor capabilities and guarantee that the piezoelectric material behaves linearly, a uniform pressure load of  0 = 1000 N/m 2 is applied to the top surface.Equation (30) is utilized to calculate the static responses of the bimorph sensor.The numerical results for the deflection and the electric potential are given with the following dimensionless units: where  11 is the stiffness constant of elastic coefficient matrix shown in (32) and the amplification factor  0 is taken as  0 = 10 10 V/m.For the purpose of comparison, a coupled 3D analysis is carried out using 20-noded hexahedral 3D piezoelectric elements (C3D20RE) with a mesh size of 40 × 20 × 10 in ABAQUS and the results from the coupled 3D FE analysis are taken as accurate.
The through-the-thickness variations of both the deflection  and the electric potential Φ at the centre of the bimorph are collected in Figure 3 for the slenderness ratio  = 5.It can be observed from Figure 3(a) that the deflection  estimated by the present method adopting different number of sublayers is constant through the thickness and it is a good approximation of the nonlinear distribution described by the coupled 3D FE analysis.The present model based on FSDT which assumes uniform deflection across the thickness cannot predict the nonlinear variation of  through the thickness.Moreover, it is noticed that by using more sublayers the deflections tend to be smaller.According to (25), a higher stiffness would be obtained with more constraints on the electric potentials across the thickness.The electric potentials induced by the plate deformation of the bimorph through the direct piezoelectric effects are shown in Figure 3(b).It  is observed that the distribution of the electric potential Φ across the bimorph thickness predicted by the present model with more than 2 sublayers is in good agreement with the nonlinear distribution predicted by the coupled 3D FE analysis.Furthermore, it is expected that with more sublayers adopted the nonlinear distribution of the electric potential Φ across the bimorph thickness can be accurately approached without introducing any higher-order electric potential assumptions.However, the conventional linear through-the-thickness electric potential model [28] would be inaccurate to predict the local electric potential response of the piezoelectric bimorph for the case of sensors.
Comparisons of the numerical results from the present model with those from the 3D FE analysis are presented in Table 1 for three typical slenderness ratios.The maximum estimating errors for the deflection are 7.29% for /ℎ = 5, 5.19 for /ℎ = 10, and 2.61% for /ℎ = 50.With 10 sublayers adopted, the discrepancy for the maximum values of the electric potential at the plate center is only 0.89% for the thin plate and 4.20% for the thick plate.Consequently, it can be concluded that the present model can provide a good approximation to the global responses such as the deflections.Furthermore, it should be noted that the local responses such as the induced electric potentials of the piezoelectric bimorph sensor can also be well predicted by the present model.

Natural Frequencies of the Piezoelectric Bimorph.
We then propose the prediction of natural frequencies of the piezoelectric bimorph plate for closed circuit condition on the top and bottom surfaces of the structure for the typical slenderness ratio  = 10.The order of Legendre polynomial is chosen as  = 4 and the mesh in Figure 4(a) is adopted.The natural frequencies of the bimorph plate predicted by the present model are shown in Table 2 in comparison to the results provided by the 3D FE analysis.A rether good agreement between the present model and the 3D FE analysis is observed (maximum error is 5.70%), which indicates that the distribution of the induced electric potential across  5.It is observed that modes 1, 3, and 6 are flexural modes, modes 2 and 5 are torsional modes, and mode 4 is an in-plane shear mode.Therefore, it should be noted that the vibration characteristics of the bimorph plate can be well addressed by the present model.The transversally isotropic piezoelectric material cannot affect the in-plane shear mode, so that the natural frequency of the fourth mode remains the same when using different number of sublayers for the PZT layer, as can be seen in Table 2.
Finally, the convergence of the present SE model is investigated.The bimorph plate is meshed into 2, 3, or 8 quadrilateral elements, as shown in Figures 4(b), 4(c), and 2, respectively.For this case, 2 sublayers are adopted.The first 2 natural frequencies are shown in Table 3 for different order of Legendre polynomials.It can be observed that the monotonic convergence rate of the present method is very fast with respect to the order of the Legendre polynomial and the algorithm is robust with respect to mesh size.

Conclusions
An efficient SE model based on the combination of an ESL approach for the mechanical displacement and a sublayer approximation for the electric potential is presented for the static and dynamic analysis of a piezoelectric bimorph.2D GLL shape functions are used to discretize the displacements.It is observed that the monotonic convergence rate of the developed SE model is very fast with respect to the order of Legendre polynomials.The capability of the present model for prediction of global and local responses such as mechanical displacements, natural frequencies, and the electric potentials across the thickness for static and dynamic processes has been verified by good agreement in numerical solutions with a 3D FE model.The advantage of the present model is that, with appropriate number of sublayers adopted, the nonlinear distribution of the electric potential across the bimorph thickness can be accurately predicted even for rather thick bimorph sensors without introducing any higher-order electric potential assumptions.It is shown that the number of sublayers adopted in the present model will not contribute to the accuracy of the global natural frequencies.Consequently, it can be further concluded that the distribution of the induced electric potential across the thickness of the bimorph sensor does not affect the global natural frequencies much.

Figure 2 :
Figure 2: Discretization of a plate and an example of spectral element.

Figure 5 :
Figure 5: Mode shapes of the first six modes of the bimorph plate for  = 10.

Table 1 :
Static responses of the piezoelectric bimorph sensor.

Table 2 :
Natural frequencies for the bimorph plate in close circuit (Hz).

Table 3 :
Natural frequencies for different order of Legendre polynomials (Hz).It can then be safely concluded that the number of sublayers adopted in the present model will not contribute to the accuracy of the global natural frequencies.By using the present SE model with 2 sublayers, the first six modes shapes of the bimorph plate are shown in Figure