Modeling of the through-the-thickness electric potentials of a piezoelectric bimorph using the spectral element method.

An efficient spectral element (SE) with electric potential degrees of freedom (DOF) is proposed to investigate the static electromechanical responses of a piezoelectric bimorph for its actuator and sensor functions. A sublayer model based on the piecewise linear approximation for the electric potential is used to describe the nonlinear distribution of electric potential through the thickness of the piezoelectric layers. An equivalent single layer (ESL) model based on first-order shear deformation theory (FSDT) is used to describe the displacement field. The Legendre orthogonal polynomials of order 5 are used in the element interpolation functions. The validity and the capability of the present SE model for investigation of global and local responses of the piezoelectric bimorph are confirmed by comparing the present solutions with those obtained from coupled 3-D finite element (FE) analysis. It is shown that, without introducing any higher-order electric potential assumptions, the current method can accurately describe the distribution of the electric potential across the thickness even for a rather thick bimorph. It is revealed that the effect of electric potential is significant when the bimorph is used as sensor while the effect is insignificant when the bimorph is used as actuator, and therefore, the present study may provide a better understanding of the nonlinear induced electric potential for bimorph sensor and actuator.

Introduction

Piezoelectric materials generate electric potentials in response to mechanical stresses, and conversely, produce mechanical movements in response to electric potentials. Therefore, piezoelectric materials can be used both as actuators and sensors, they transform electrical energy into mechanical energy, and vice versa. To achieve practically meaningful actuation and sensing capabilities, a piezoelectric bimorph consisting of two piezoelectric layers is widely used [1,2]. A broad range of electromechanical applications have been reported, such as electroacoustic transducers [3,4], medical devices [5], microcantilever biosensors [6], and atomic force microscope (AFM) cantilevers [7]. However, before piezoelectric bimorphs can be utilized in all these applications, it is first necessary to investigate both the global responses and the local responses, e.g., the deflection and the distribution of the electric potential across the thickness.

There have been many theories and models developed for analyzing piezoelectric bimorph structures with emphasis on approximating the mechanical displacement and electric potential. By carrying out exact 3-D analytical solutions for the simply supported piezoelectric plate [8,9], it is shown that the distribution of the electric potential across the thickness is nearly quadratic. This implies that the assumption of linear distribution of the electric potential across the thickness adopted by many numerical models [10,11] cannot address this nonlinear electric potential. Since exact 3-D analytical solutions are not available for more general cases of loading and boundary conditions, the introduction of the finite element (FE) method is desirable. A considerable amount of literature has been published on the FE analysis of piezoelectric smart structures [12-14]. Among these works, the simplest and often used model is the equivalent single layer (ESL) model in which the displacement and strain functions are assumed to be continuous through the thickness. There are two main kinds of theories used for ESL models. One is the classical laminated plate theory (CLPT) [15,16], and the other one is the shear deformation theory, which branches out into first-order shear deformation theory (FSDT) [17,18] and higher order shear deformation theory (HSDT) [19,20]. 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 the thickness. To overcome this shortcoming, the FE model using the layer-wise theory [21-24] or the sublayer theory [2,25-28] has been recommended. In the latter case, the piezoelectric layer is divided into appropriate number of thin sublayers. For each of these sublayers, a linear electric potential distribution across the plate thickness is assumed. It is further expected that the quadratic distribution of the electric potential across the plate thickness can be accurately approached with more sublayers adopted.

Generally, accurately simulation of the local responses of the piezoelectric bimorph structures would inevitably lead to a very dense FE mesh when using the FE method. Hence, conventional FE simulation becomes computationally very inefficient. A more efficient method is the spectral element (SE) method which combines the geometric flexibility of FE method with the high accuracy of the pseudo spectral method. This method was first presented by Patera in the mid 1980s [29]. 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 Cheybysev polynomials, in the shape functions. The SE method results naturally in diagonal mass matrices which is a distinct advantage over traditional FE method especially for transient analysis. Moreover, to have an accurate simulation with the conventional FE method, a mesh with a large number of elements and degrees of freedom (DOFs) 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. The SE method has been widely applied to many engineering problems related to acoustics, fluid dynamics and seismology [30-35]. Recently, the SE method has been extensively used to investigate the wave propagation problems for the purpose of damage detection in structures [36,37]. However, according to the authors' best knowledge, the SE method has not been previously used for accurately modeling of the through-the-thickness electric potentials for piezoelectric bimorphs.

For the purpose of accurately representing the mechanical displacement and the electric potential, a reasonable choice is to use the ESL model for the mechanical variables and the layer-wise theory or the sublayer theory for the electric variables. In the present work, we attempt to combine the merits of the SE method and the sublayer model. More specifically, the mechanical variables, i.e., the displacements, are described based on FSDT. The electrical variables, i.e., the potentials, are described using the sublayer model. SE method is then utilized to deduce the governing equations. Legendre orthogonal polynomials are adopted in the interpolation function to improve the accuracy. To validate the effectiveness and the capability of the present model, numerical simulations for a simply supported piezoelectric bimorph with two different load cases, i.e., a uniform pressure load applied to the top surface and a uniform potential applied to the top and bottom surfaces, are carried out. The results obtained by the present approach are then compared to those coming from the coupled 3-D FE simulations using ABAQUS. The comparisons show the good accuracy and efficiency of SE method for modeling of the through-the-thickness electric potentials of the piezoelectric bimorph.

Mathematical Formulation

Constitutive Relationships, Displacement and Strain

A piezoelectric bimorph made of two identical PZT-4 piezoelectric layers, which has been investigated by Fernandes [1], is considered here. The PZT-4 layer is assumed to behave in a linear orthotropic manner with small displacements and strains. As depicted in Figure 1, both piezoelectric layers have the same thickness 0.5 h and are poled in the same direction. The x-y plane of the coordinate system x-y-z coincides with the middle plane of the bimorph, and the z axis is defined normal to the middle plane following the right-hand rule. This work aims to investigate the problem of a simply supported piezoelectric bimorph under a uniform pressure load or an applied electric potential in the framework of linear theory of piezoelectricity. Assuming the PZT-4 layers work under isothermal conditions, the pyroelectric effects and thermomechanical couplings are not taken into account. Consequently, a linear constitutive relationship addressing both the direct and converse piezoelectric effects is utilized for the analysis of the piezoelectric bimorph, which can be written as: where σ = [σ]T and ε = [ε]T represent stress vector and strain vector, respectively. E = [E]T, the electric field vector, D = [D]T, the electric displacement vector, c, the elastic coefficient matrix, g, the dielectric coefficient matrix, and e, the piezoelectric stress coefficient matrix.

Figure 1.

Geometry of a piezoelectric bimorph.

An ESL model adopting the FSDT is adopted to describe the mechanical displacement. The displacement field of a piezoelectric bimorph based on FSDT takes on the form [17,18]: where u̅, υ̅, w̅ denote the displacements of an arbitrary point on the mid plane z = 0, α̅ and −β̅ denote the rotations of a transverse normal about the y and x axes, respectively. In the FSDT, the transverse shear strains are assumed to be constant with respect to the thickness coordinate. The constant state of transverse shear strains across the thickness is a gross approximation of the true strain field, which is at least quadratic through the thickness.

We define: where U is the displacement vector, and U̅ is a generalized displacement vector. Then Equation (2) can be written in matrix form as: where:

The infinitesimal strain components associated with the displacements are given by: where L is the derivation operator defined as:

Approximations for Displacements

The Legendre polynomials based SE method can be described as follows: the bimorph is firstly discretized using a set of non-overlapping rectangular elements, as in the traditional FE method. Each rectangular element, denoted by Ωe, is then mapped to a reference element, denoted by Ωref: ξ ∈ [−1,1]× η ∈ [−1,1], using an invertible local mapping. The discretization procedure is illustrated in Figure 2.

Figure 2.

Discretization of a plate and an example of spectral element.

Subsequently, a set of nodes, denoted by ξ, η are defined in the local coordinate system ξ−η of the reference element Ωref as roots of the following polynomial expression: where P is the N-th order Legendre polynomial. In fact, the nodes are the 2-D Gauss-Lobatto-Legendre (GLL) points. In contrast to the classical FE method, the distribution of nodes is irregular, as shown in Figure 2. In the current formulation, the 5-th order Legendre polynomial is chosen, hence 36 nodes can be specified in the reference element Ωref, as depicted in Figure 2.

The 1-D shape functions at the 1-D GLL points ξ are defined as [36]: An important property of these interpolation functions is the discrete orthogonality expressed as: where δ denotes the Kronecker delta. The 2-D shape functions are constructed as a tensor product of the 1-D ones:

Figure 3 shows two examples of the 2-D shape functions which indicate that each shape function has the value 1 at one node and vanish at all other nodes.

Figure 3.

Selected shape functions for a 36-node spectral element. (a) Ψ32(ξ,η); (b) Ψ45(ξ,η).

Coordinates x and y within each Ωe may be uniquely related to ξ and η upon the invertible mapping: where x and y denote the coordinates of x and y, respectively, of the element nodes ξ, η. The generalized displacements u̅, υ̅, w̅, α̅ and β̅ over an reference element Ωref are discretized by the 2-D shape functions as: where u̅, υ̅, w̅, α̅ 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 Equations (14) into Equation (15) yields: where N is the displacement shape function matrix which can be expressed as: with: where I5×5 is a 5×5 identity matrix. Substituting Equation (17) into Equation (7) yields: where B is strain-displacement matrix which can be written as:

Approximations for Electric Potential

For the purpose of accurately modeling the distribution of the electric potential across thickness, each layer of the piezoelectric bimorph is subdivided mathematically into n thinner sublayers. As shown in Figure 4, the sublayers are numbered in top-to-bottom order. The z coordinates of the top and bottom surfaces of the i-th sublayer are denoted by z and z−1, respectively. In each sublayer, the distribution of the electric potential ϕ(z) is assumed to be linear across the thickness such that: where is the interpolation function and Φ̃ is a column matrix composed of the electric potentials at the top and the bottom surfaces of the i–th sublayer, which can be expressed as:

Figure 4.

A sublayer model for a piezoelectric bimorph.

In this way, the assumption of linear distribution of electric potential across the thickness is used not in the whole piezoelectric layer, but in each sublayer instead. As a result, the electric potential is approximated as piecewise linear across the thickness and it is expected that the quadratic distribution of the electric potential across the bimorph thickness can be approached with more sublayers adopted. As mentioned before, the mechanical displacement field is approximated using ESL model based on FSDT and the piezoelectric bimorph is discretized using 2-D mesh. To keep the compatibility, each sublayer of the piezoelectric layer is also discretized using the same mesh. Consequently, an element potential vector Φe is then introduced in the spectral plate finite element, which is defined as:

The surface potential of the sublayer, ϕ, is assumed to be constant over the element and ϕ0, ϕ1, ⋯ ϕ2 are treated as elemental 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 x, y. Thus the present method combines an ESL theory for the displacement 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(z) is the electric field of the i-th sublayer, is the electric field-potential matrix, given by:

Governing Equations

By applying Hamilton's principle, the elementary dynamic equations for the piezoelectric bimorph plate can be obtained: where denotes the element mass matrix; , mechanical stiffness matrix; and the piezoelectric coupling matrices; the dielectric permittivity matrix; the vector of externally applied force; and the vector of externally applied charge, respectively: where ρ is the mass density, P is the surface force vector, q is the surface charge density vector. J is the well-known Jocobian matrix of the mapping (13) which is defined as:

The GLL integration rule is then used to calculate the characteristic matrices and the nodal force vector in Equation (28) at the elemental level [36]. In this study, the interface between the two PZT layers is grounded. Two sets of electric boundary conditions are considered, i.e., (1) sensor function with the top and bottom surfaces grounded and a uniform pressure load of S0 = 1,000N/m2 applied to the upper surface, and Equation (2) actuator function with an electric potential of V0 = 50V applied to the top and bottom surface of the bimorph. By applying the electric boundary conditions, the DOFs for the electric potential are condensed out such that Equation (28) is finally of the form: where denotes the mechanical stiffness matrix induced by the electromechanical coupling of PZT-4 layer, and denotes 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 , , , and ; q is the global nodal displacement vector. Since the DOFs for the sublayer electric potentials have been condensed out, this approach will not result in a large number of potential DOFs. For the purpose of static analysis, the governing equations in Equation (37) reduces to:

Numerical Results

In this section, the derived SE model is converted into a numerical code and case studies are carried out to validate the effectiveness and the capability of the present model for predicting both the global responses and the local responses, i.e., the deflections of the bimorph and the distribution of the electric potential across the bimorph thickness. A simply supported rectangular piezoelectric bimorph shown in Figure 1, which has been investigated by Fernandes [1], is considered here. The material constants of PZT-4 are given as:

The length a and width b of the bimorph are 25 mm and 12.5 mm respectively. Two values of slenderness ratio, S = a/h = 5 and, S = 50 which represent the 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 3 is used in this work. Two load cases corresponding respectively to sensor function and actuator function are considered. To overcome the ill condition problem resulted from the huge difference of the element values of and in magnitudes, Equation (28) is rewritten using dimensionless variables. Consequently, the numerical results for the deflection and the electric potential are given in dimensionless units as: where the amplification factor E0 is taken as E0 = 1010V/m. For the purpose of comparison, a coupled 3-D analysis is carried out using 20-noded hexahedral 3-D piezoelectric elements (C3D20RE) with a mesh size of 40×20×10 in ABAQUS and the results from the coupled 3-D FE analysis are taken as accurate.

Sensor Function

For this case a uniform pressure load of S0 = 1,000N/m2 is applied to the upper surface and the bimorph is used as a sensor with the top and bottom surfaces grounded. The variations of both the deflection W and the electric potential Φ across the bimorph thickness at the centre of the bimorph plate (x = 0.5a, y = 0.5b) for the slenderness ratio S = 5 and S = 50 are shown in Figures 5 and 6, respectively. It can be observed from Figure 5a and Figure 6a that the deflection W 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 3-D analysis. The present model based on FSDT with assumption of uniform deflection through the thickness cannot predict the nonlinear variation of W through the thickness. The electric potentials induced by the deformation of the bimorph through the direct piezoelectric effects are shown in Figures 5b and 6b. It is observed that the distribution of the electric potential Φ across the thickness provided by the present approach with more than 2 sublayers is in good agreement with the nonlinear distribution predicted by the coupled 3-D analysis. Furthermore, it is expected that with more sublayers adopted the quadratic 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 electric potential assumption [38] will result in an inaccurate prediction of the local electric potential response for the case of sensors. The curves in Figures 5 and 6 are symmetrical with respect to the interface between the two PZT layers. It should be highlighted that although the present method cannot predict accurately the distribution of W across the bimorph thickness, it may be able to provide good approximate results for Φ with appropriate number of sublayers for both thick and thin bimorph plates.

Figure 5.

Bimorph sensor of S = 5 under pressure load. (a) Dimensionless deflection; (b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 5 (triangles) and present model with n = 5 (small circles).

Figure 6.

Bimorph sensor of S = 50 under pressure load. (a) Dimensionless deflection; (b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 2 (triangles) and present model with n = 10 (small circles).

Actuator Function

To achieve practically meaningful actuation capabilities and guarantee that the piezoelectric material behaves linearly, an electric potential of V0 = 50V is applied to the top and bottom surfaces of the bimorph with intermediate electrode grounded. The through-the-thickness variations of W and Φ at the centre of the plate for S = 5 and S = 50 are shown in Figures 7 and 8, respectively.

Figure 7.

Bimorph actuator of S = 5 under potential load. (a) Dimensionless deflection; (b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 2 (triangles) and present model with n = 10 (small circles).

Figure 8.

Bimorph actuator of S = 50 under potential load. (a) Dimensionless deflection; (b) Dimensionless electric potential. 3-D FE analysis (full line), present model with n = 2 (triangles) and present model with n = 10 (small circles).

Once again, the present model based on FSDT cannot predict accurately the through-the-thickness distribution of W. Similar to the previous observation, the constant deflection W through the thickness calculated by the present method adopting different number of sublayers is a good approximation of the nonlinear distribution provided by the coupled 3-D analysis. It is noticed that as the sublayer number increases a smaller deflection is obtained which is also pointed out by Wang [2]. The electric potentials at the centre of the plate are plotted in Figure 7b and Figure 8b for the slenderness ratio S = 5 and S = 50, respectively. It can be observed that the almost linear distribution of Φ across the thickness predicted by the present method for both thick and thin bimorph plates is in excellent agreement with the coupled 3-D analysis, indicating that the nonlinear induced electric potential is insignificant compared to the externally applied potential. Consequently, the conventional linear electric potential assumption [38] may be accurate enough to calculate the local electric potential response for the case of actuators.

Conclusions

The present work aims to develop an efficient SE model with electric potential DOFs for the static electromechanical response of a piezoelectric bimorph. The approach is the combination of an ESL model based on FSDT for the mechanical displacement with a sublayer model based on the piecewise linear approximation for the electric potential. 2-D GLL shape functions are used to discretize the displacements and then the governing equation of motion is derived following the standard SEM procedure. By applying the electric boundary conditions, the DOFs for the electric potential are condensed out such that the present model will not result in a large number of potential DOFs.

Numerical simulations based on the present model are carried out for two different load cases, i.e., a uniform pressure load applied to the top surface and a uniform potential applied to the top and bottom surfaces. To validate the effectiveness and the capability of the present model for investigation of both global and local response of the piezoelectric bimorph, the numerical results thus obtained are compared to those from 3-D analysis using ABAQUS. The results indicate that the deflection W estimated by the present method is a good approximation of the nonlinear distribution predicted by the coupled 3-D analysis. It is further shown that the present model provides very accurate prediction for the electric potential distributions across the bimorph thickness even for rather thick bimorph plate without introducing any higher-order electric potential assumptions. It is also revealed that the conventional linear electric potential model is accurate enough to predict the local electric potential response for the case of actuators. This observation consists with the previous findings proposed by Yang [39]. One of the limitations is that the deflection W across the thickness is constant. Nevertheless, it is accurate enough to investigate the global response of the piezoelectric bimorph. The present work is important for researchers to better understand the nonlinear induced electric potential for bimorph sensor and actuator. An important extension of the present research is to study the vibration characteristics of the piezoelectric bimorph based on SE method. The influence of the induced stiffness matrix on the natural frequencies of the bimorph plate under various electric boundary conditions is to be investigated. The convergence study of the present model with respect to the order of the Legendre polynomial is also a practical and interesting problem to be conducted.