An electromechanical Reissner–Mindlin model for laminated piezoelectric plates

Composite Structures 88 (2009) 394–402

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

An electromechanical Reissner–Mindlin model for laminated piezoelectric plates Lin Liao, Wenbin Yu * Mechanical and Aerospace Engineering Department, Utah State University, Logan, 84322-4130, USA

a r t i c l e

i n f o

Article history: Available online 1 May 2008 Keywords: Variational asymptotic method Piezoelectric plate Reissner–Mindlin

a b s t r a c t An electromechanical Reissner–Mindlin model is constructed for laminated piezoelectric plates using the variational asymptotic method. This model is applicable to laminates without prescribed electric potential through the thickness. Taking advantage of the smallness of the plate thickness, we rigorously split the original 3D piezoelectricity problem into a 1D through-the-thickness analysis and a 2D plate analysis, and both are fully-coupled electromechanical analyses. Examples of single layer and multi-layer plates have been used to demonstrate the accuracy and application of this model. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Considerable amount of research and investigation have been done to study laminated piezoelectric plates. The exact solutions of laminated piezoelectric plates can be obtained analytically for a few cases of ideal material type, geometry, and boundary conditions [1–3]. For general cases, we can use the three-dimensional (3D) multiphysics finite element method (FEM) to find numerical solutions [4]. However, it is labor intensive to prepare the 3D multiphysics FEM model for a laminated piezoelectric plate, especially if it is composed of many anisotropic layers with different orientations. Moreover, the prohibitive computational cost of 3D FEM can only be justified for detailed analyses and prevents its use in design and overall simulation of realistic engineering systems involving such components. Numerous two-dimensional (2D) plate models have been constructed in order to simplify the analysis of piezoelectric plates, which generally start from some assumptions for the through-the-thickness distribution of the 3D field quantities. According to the way the assumptions applied, piezoelectric plate models can be classified as equivalent single-layer models [5–17] if the assumptions are applied to the entire structure and layerwise models [18–21] if the assumptions are applied to each layer. In the past investigation, single-layer theory of mechanical fields is usually combined with layerwise approximation of electric potential [22–25], and 2D finite element is incorporated in developing these plate theories [22,23,26]. In addition, considerable efforts have been dedicated to the study of vibration suppression, shape control, and buckling enhancement, and optimization of composite plates with embedded or surface bonded piezoelectric actuators/ sensors [27–29]. These models have two main disadvantages: * Corresponding author. Tel.: +1 435 7978246; fax: +1 435 7972417. E-mail address: [email protected] (W. Yu). 0263-8223/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2008.04.017

(1) the a priori assumptions which are natural extensions from isotropic, homogeneous structures cannot be easily justified for highly heterogeneous and anisotropic structures such as laminated piezoelectric plates; (2) there is no rational way for the analyst to determine the loss of accuracy and which refinement (i.e. singlelayer versus layerwise, first-order versus higher-order) should be used for a reasonable tradeoff between accuracy and efficiency. Extensive researches have been done on piezoelectric plates with electroded face surfaces and interfaces between layers [30] so that the electric potential can be prescribed at a certain point through the thickness. However, it is also possible that no electric potential is prescribed at any point through the thickness, say the electrodes are coated on the lateral boundary of the plate. The research works on this electrode arrangement are associated with transducers and resonators, which oftentimes belong to the area of electrical engineering or physics. It is noted that the lateral field excitation can also be produced by placing two electrodes on the same face surface of the plate [31–34], which might be one of the reasons that study on the case of electroded lateral boundary has not received much attention. Most of the investigations for this type of loading have been devoted to the study on vibration modes, depolarizing-field effect, electromechanical coupling coefficients, and admittance [35–38]. Published works related to the mechanics of piezoelectric composite plates with electroded lateral edges are rare. In particular, static analysis of displacements, strains, and stresses has not been found in the literature. Even for prescribed electric load on the lateral boundary, the electric field is often specified as a constant instead of a field generated by prescribed electric potential along the electroded edges [39]. The focus of this paper is to construct a generalized Reissner– Mindlin model with both electric and mechanical measures for laminated piezoelectric plates without prescribed electric potential through the thickness. The variational asymptotic method (VAM)

395

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402

[40] will be applied to mathematically split the coupled 3D electromechanical analysis into a series of one-dimensional (1D) electromechanical through-the-thickness analyses and 2D electromechanical analyses over the reference plane. The coupled through-the-thickness analyses calculate both the mechanical warping and the electric warping, the meaning of which will be defined later. The results of the coupled through-the-thickness analyses are generalized 2D electric enthalpy, which can be transformed into the form of a generalized Reissner–Mindlin model. The coupled through-the-thickness analyses can also provide a set of recovery relations to reproduce 3D distribution of the electromechanical fields within the structure based on the results obtained from the 2D electromechanical analysis. One of the major applications of 2D plate models developed is the computation of resonant frequencies of piezoelectric devices [41]. 2. Three-dimensional formulation The Hamilton’s extended principle for a piezoelectric composite plate can be written as

Z

t2

½dðK  HÞ þ dWdt ¼ 0

ð1Þ

t1

where t1 and t 2 are arbitrary fixed times, K and H are the kinetic energy and electric enthalpy, respectively, dW is the virtual work of applied loads and electric charges (if exist). The bar is used to indicate that the virtual work need not be the variations of functionals. For piezoelectrics, the electric enthalpy can be expressed as

H¼

1 2

Z

ð CT : C E : C  2E  e : C  ET  eC  EÞdV

ð2Þ

V

where C E is the elastic tensor at constant electric field, C is the strain tensor, e is the piezoelectric tensor, E is the electric field vector, eC is the dielectric tensor at constant strain field, and V is the space occupied by the structure. Although we are focusing on plates made of piezoelectrics, the present formulation is equally applicable to smart structures made of other active materials characterized by a constitutive model with the same mathematical structure. As sketched in Fig. 1, a point in the plate can be described by its Cartesian coordinates xi , where xa are two orthogonal lines in the reference plane and x3 is the normal coordinate. (Here and throughout the paper, Greek indices assume values 1 and 2 while Latin indices assume 1, 2, and 3. Repeated indices are summed over their range except where explicitly indicated.) Letting bi denote the

B b

^rðx1 ; x2 ; x3 Þ ¼ rðx1 ; x2 Þ þ x3 b3

ð3Þ

where r is the position vector from O to the point located by xa on the reference plane. When the reference plane of the undeformed plate coincides with its middle plane, we have

h^rðx1 ; x2 ; x3 Þi ¼ hrðx1 ; x2 Þ

ð4Þ

where the angle-brackets denote the definite integral through the thickness of the plate, denoted as h. The position vector of any ^ can be represented as material point in the deformed plate R

^ 1 ; x2 ; x3 Þ ¼ Rðx1 ; x2 Þ þ x3 B3 ðx1 ; x2 Þ þ wi ðx1 ; x2 ; x3 ÞBi ðx1 ; x2 Þ Rðx

ð5Þ

where R is the position vector of the reference plane for the deformed plate, Bi is an orthonormal triad for the deformed configuration, and wi are warping functions introduced to accommodate all possible deformation other than those described by R and Bi . To en^ in terms of R, Bi , and wi , we need to sure a unique expression for R introduce six constraints. We can define R to be the average position through the plate thickness, which implies that the warping functions must satisfy the following three constraints

hwi ðx1 ; x2 ; x3 Þi ¼ 0

ð6Þ

Another two constraints can be specified by choosing B3 as the normal to the reference plane of the deformed plate. It should be noted that this choice has nothing to do with the famous Kirchhoff hypothesis, in which no local deformation of the transverse normal is allowed. Here, we have accommodated all possible deformation using the warping functions. Because Ba can freely rotate around B3 , we can introduce the last constraint as

B1  R;2 ¼ B2  R;1

ð7Þ

Based on the concept of decomposition of rotation tensor [42], we can obtain 3D strains valid for small local rotations using

1 2

Cij ¼ ðF ij þ F ji Þ  dij

ð8Þ

where dij is the Kronecker symbol, and F ij is the mixed-basis component of the deformation gradient tensor such that

F ij ¼ Bi  Gk gk  bj

ð9Þ

with Gk as the covariant base vectors of the deformed configuration and gk as the contravariant base vectors of the undeformed configuration. The 2D generalized strains eab and K ab can be defined as

R;a ¼ Ba þ eab Bb

ð10Þ

b

Bi;a ¼ ðK ab Bb  B3 þ K a3 B3 Þ  Bi

r

Using Eq. (8) along with Eqs. (9), (10), (5) and (3), we can express the 3D strain field Cij in terms of eab , K ab , and wi . For geometrically nonlinear analysis, we can assume that the strains are small compared to unity and warpings are of the order of strain or smaller. Neglecting higher-order terms, one can express the 3D strain field as

b B

B

u

r

unit vector along xi for the undeformed plate, one can then describe the position of any material point in the undeformed configuration by its position vector ^r from a fixed point O, such that

1

2

R

B

3

1

2

Ce ¼ e þ x3 j þ I1 wk;1 þ I2 wk;2 2Cs ¼ w0k þ e1 w3;1 þ e2 w3;2

3

Ct ¼

w03

R oðÞ where ðÞ0 ¼ ox , ðÞk ¼ bðÞ1 ðÞ2 cT 3

Ce ¼ bC11 2C12 C22 cT 2Cs ¼ b2C13 2C23 cT Ct ¼ C33 Fig. 1. Schematic of plate deformation.

e ¼ be11 2e12 e22 cT

j ¼ bK 11 K 12 þ K 21 K 22 cT

ð11Þ

396

and

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402

the top and bottom surfaces, respectively, Q ¼ Q i Bi is the applied b is the applied electric tractions along the lateral edges, and D charges along the lateral surfaces. ðÞþ or ðÞ denote the correspondb ing quantity evaluated at the top or bottom surface, respectively. d R is the Lagrangian variation of the displacement field, such that

2

3 2 3 1 0 0 0     1 0 6 7 6 7 I1 ¼ 4 0 1 5 I2 ¼ 4 1 0 5 e1 ¼ e2 ¼ 0 1 0 0 0 1

Up to this point, we have obtained an exact formulation for the kinematics of the mechanical field of a piezoelectric plate, which is the same as conventional composite plates [43]. However, a complete description of the piezoelectric plate requires not only the mechanical field but also the electric field, which is characterized by the electric potential, /ðxi Þ. Following the reasoning of [44], if there is no prescribed potential at any point through the thickness, the electric potential can be expressed as using the following change of variables:

/ðx1 ; x2 ; x3 Þ ¼ Uðx1 ; x2 Þ þ gðx1 ; x2 ; x3 Þ

ð12Þ

where Uðx1 ; x2 Þ is defined as the average of /ðx1 ; x2 ; x3 Þ through the thickness, which implies

hgðx1 ; x2 ; x3 Þi ¼ 0

ð13Þ

Here we term U as the 2D electric potential, g as the electric warping function which is asymptotically smaller than U. The 3D electric field is defined as the negative gradient of the 3D electric potential, which is

Es ¼ E2D  g;a

Et ¼ g 0

ð14Þ

with E2D ¼ U;a as the 2D electric field. To obtain the kinetic energy, we can calculate the absolute velocity of a generic point in the structure by taking a time derivative of Eq. (5), such that

~ ðn þ wÞ þ w _ v¼V þx

ð15Þ

_ is the partial derivative with respect to time, V is the absowhere ðÞ lute velocity of a point in the deformed reference plane, x is the e forms an inertial angular velocity of Bi bases, and the notation ðÞ e ij ¼ eijk ðÞ using antisymmetric matrix from a vector according to ðÞ k the permutation symbol eijk . In Eq. (15), the symbols v; V; x; w denote column matrices containing the components of corresponding vectors in Bi bases, and n ¼ b0 0 x3 cT . The kinetic energy of a plate can be obtained by

1 K¼ 2

Z

T



qv vdV ¼ K2D þ K

ð16Þ

ð21Þ

where the virtual displacement and rotation are defined as

dqi ¼ dR  Bi

dBi ¼ ðdw2 B1 þ dw1 B2 þ dw3 B3 Þ  Bi

ð22Þ

where dqi and dwi are components of the virtual displacement and rotation in the Bi system, respectively. Since the warping functions are small, one may safely ignore products of the warping and virtual b and obtain the virtual work due to applied loads as rotation in d R

dW ¼ dW2D þ dW

ð23Þ

where

dW2D ¼

Z

ðfi dqi þ ma dwa  ðDþ þ D ÞdUÞdX Z b UÞds þ ðhQ i idqi þ hx3 Q a idwa  h Did X

dW ¼

Z

ð24Þ

oX

ðhPi dwi i þ si dwþi þ bi dwi  Dþ dg þ  D dg  ÞdX Z b þ hQ i dwi  Ddgids

X

ð25Þ

oX

with the generalized forces fi and moments ma defined as

fi ¼ hP i i þ si þ bi

h ma ¼ hx3 Pa i þ ðsa  ba Þ 2

ð26Þ

The second integration in Eq. (25) accounts for virtual work done along the lateral boundaries of the plate due to the mechanical warping and the electric warping. This term is necessary for the edge-zone problem, which is an important subject in its own right and beyond the scope of the present paper. Therefore, we will drop this term hereafter. With the knowledge of Eqs. (16) and (23), the Hamilton’s extended principle in Eq. (1) becomes

Z

t2

½dðK2D þ K  HÞ þ dW2D þ dW dt ¼ 0

ð27Þ

t1

V

where q is the mass density and

Z

1 fnV þ xT ixÞdX  V T V þ 2xT l ðl ð17Þ 2 X Z h i 1 e w þ wÞ ~ w þ wÞ ~ nÞT ðx ~ w þ wÞ _ T ðx _ þ 2ðV þ x _ dV K ¼ q ðx 2 V ð18Þ

K2D ¼

 ; l n, and i are inertial constants commonly used in plate where l dynamics, which can be trivially obtained through simple integrals over the thickness as:

2

hx23 qi

l ¼ hqi ln ¼ b0 0 hx3 qicT i ¼ 6 4 0

0

0

0

0

0

7 hx23 qi 0 5

ð19Þ

Z

b þ s  dR b þ þ b  dR b   Dþ d/þ  D d/ ÞdX ðhP  d Ri Z b  Dd/ids b þ hQ  d R

So far, we have presented a 3D formulation for the electromechanically coupled problem of piezoelectric plates in terms of 2D displacements (represented by R  r) and rotations (represented by bi and Bi ) and 3D warping functions (wi and g). If we attempt to solve this problem directly, the difficulty associated with the original 3D model will arise. Fortunately, we can use the variational asymptotic method [40] to solve wi and g through an asymptotical analysis of the variational statement in Eq. (27) in terms of small parameter of the aspect ratio inherent in a plate structure.

3. Dimensional reduction

3

If no charge density exists inside the body, the virtual work of the smart plate can be calculated as

dW ¼

b ¼ dqi Bi þ x3 dB3 þ dwi Bi þ wj dBj dR

X

ð20Þ

oX

where oX denotes the boundary of the reference plane, P ¼ Pi Bi is the applied body force, s ¼ si Bi ; b ¼ bi Bi are tractions applied on

The dimensional reduction from the original 3D formulation to a 2D plate model can only be performed approximately. One way to do it rigorously is to take advantages of the small parameters of plates, h=l  1, with l as the characteristic wavelength of the deformation of the reference plane. The strain and electric field are also small if we are only interested in a geometrically nonlinear but physically linear 2D theory, i.e.

e  hj  g  v=lEi

ð28Þ

l is characteristic magnitude of the elastic constants, v is the characteristic magnitude of the piezoelectric coefficients, and g is the order of the maximum strain in the plate. From the plate equations

397

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402

of equilibrium, we can estimate the orders of the following quantities corresponding to the order of strains:

hP3  s3  b3  lðh=lÞ2 g hPa  sa  ba  lðh=lÞg

ð29Þ

Dþ  D  vðh=lÞ2 g

where k3 , kk , and k4 are the corresponding Lagrange multipliers to enforce the constraints for w3 ; wk , and g, respectively. The natural boundary conditions on the face surfaces and natural continuity conditions on the interfaces enable us to simplify the Euler– Lagrange equation in Eq. (35) to be: T 0T T 0T ðe þ x3 jÞT C et þ w0T k C st þ w3 C t  E2D ec þ g et ¼ 0 T 0T T T 0T T ðe þ x3 jÞT C es þ w0T k C s þ w3 C st  E2D es þ g ea ¼ 0

3.1. Analytical dimensional reduction

T

As the first step, a classical model approximating the original energy up to the zeroth order – lg2 will be constructed. If we choose the characteristic scale of change of the displacements and warping functions in time in such a way that K is much smaller than other terms in Eq. (27), which is valid for most general applications [44], then according to VAM, we have the following variational statement instead:

Z

t2

½dðK2D  HÞ þ dW2D þ dW dt ¼ 0

ð30Þ

t1

in which the electric enthalpy H can be rewritten in terms of elecR tric enthalpy density per unit area H, as H ¼ X HdX, in which H is expressed analytically as

9 2 8 *> Ce >T C e = < 1 6 CT H¼ 2 Cs 4 es > > 2 : ; Ct C Tet 9T 2 8 *> Ce > eb = < 6  2 Cs 4 es > > ; : Ct ec * T  ds Es 1  T 2 Et det

C es Cs

9 38 + C et > = < Ce > 7 C st 5 2Cs > > ; : Ct Ct

C Tst 3 eet  + 7 Es ea 5 Et et  + Es det Et dt

ð32Þ

t1

ð33Þ

along with the constraints in Eqs. (6) and (13), in which the electric enthalpy density associated with H0 is expressed as

   1 0T 0 0 H0 ¼ ðe þ x3 jÞT C es w0k þ C et w03 þ w0T C w þ C w w st t k 3 3 2 3 h i 1 0T 0T 0 0 þ ðe þ x3 jÞT eet þ w0T k ea þ w3 et g þ wk C s wk 2   1 0T 0T 0T T 0 E e þ w e  g d þ d g g  w0T 2D t et k s 3 c 2

 g dt ¼ 0

T  1 1 w0T þ ET2D eT 3 ¼ ðe þ x3 jÞ C et C t c Ct 1 ÞT eet dt

0

g ¼ ðe þ x3 j

þ

ð37Þ

 1 ET2D det dt

where the starred quantities are algebraic expressions in terms of the original material constants including elastic, piezoelectric, and dielectric coefficients listed in Table 1. Warping functions can be solved from Eq. (37) with the continuity conditions on the interfaces, along with the constraints in Eqs. (6) and (13). Plugging Eqs. (11) and (14) along with the solved warping functions into the original functional in Eq. (31), one will find out that the zeroth-order approximation of the electric enthalpy density asymptotically correct through the order of lg2 can be expressed in terms of e, j, and E2D as:

 T  e A

B

BT

D

j



e

j



 ET2D d2D E2D  2

 T   e e

j

ej

E2D

where A; B; D are the 2D stiffness matrices, and e ; ej are the 2D piezoelectric constants, and d2D is the 2D dielectric constants. These matrices can be evaluated using the following expressions:





 B ¼ x3 C e D ¼ x23 C e A ¼ C e D E  1   1 2e ¼ 2eb  C es C 1 s es  C es C s es  C et ec C t D E T 1 T 1 1  C  þ eet det dt þ eet det dt et ec C t D E  1   1 Þ 2ej ¼ x3 ð2eb  C es C 1 s es  C es C s es  C et ec C t D E T 1 T 1 1 þ eet det dt þ eet det dt Þ  x3 ðC  et ec C t D E T 1  T  1 d2D ¼ ds þ eTs C 1  det det dt s es þ ec ec C t

ð39Þ

where 1

T T 1 C e ¼ C e  C es C 1 þ eet eT et dt s C es  C et C et C t

ð40Þ

3.2. Numerical dimensional reduction using the finite element method

ð34Þ

Following the usual procedure of calculus of variations and considering that the warping functions are essentially different functions for each layer, we can obtain the Euler–Lagrange equations for each layer:

h i0 T 0T T 0T ðe þ x3 jÞT C et þ w0T ¼ k3 k C st þ w3 C t  E2D ec þ g et h i0 T T T 0T T 0T T ðe þ x3 jÞ C es þ w0T ¼ kk k C s þ w3 C st  E2D es þ g ea h i0 T 0T 0T ðe þ x3 jÞT eet þ w0T ¼ k4 k ea þ w3 et þ E2D det  g dt

ð36Þ

0T

ð38Þ

It can be seen that wi and g only appear in the electric enthalpy in Eq. (32) and can be solved from

dH0 ¼ 0

þ

ET2D det

T  1 T T 1 w0T k ¼ ðe þ x3 jÞ C es C s þ E2D es C s

ð31Þ

t2

½dðK2D  H0 Þ þ dW2D dt ¼ 0

þ

w0T 3 et

Eq. (36) can be used to solve for w0i and g 0 algebraically for each layer as:

2H2D ¼

where C e , C es , C et , C s , C st , C t are the corresponding partition matrices of the 6  6 elastic material matrix at constant electric field, eb , eet , es , ea , ec , et are the corresponding partition matrices of the 6  3 piezoelectric coefficient matrix, and ds , det and dt are the corresponding partition matrices of the 3  3 dielectric coefficient matrix at constant strain. To solve the zeroth-order warping, the variational statement in Eq. (30) becomes

Z

ðe þ x3 jÞ eet þ

w0T k ea

ð35Þ

Although we can obtain the zeroth-order warping analytically, the derivation becomes intractable for solving the first-order warping functions for construction of refined models. Hence, we will use Table 1 Expressions of starred quantities 

  T C  es ¼ C es þ eet ea =dt    T C  et ¼ C et þ eet et =dt    T e c ¼ ec  et det =dt    T e s ¼ es  ea det =dt   eet ¼ eet  C es C 1 e a  C et et =C t s et ¼ et  C Tst C 1 e a s ea ¼ ea  C st et =C t

C es ¼ C es  C et C Tst =C t C et ¼ C et  C es C 1 s C st ec ¼ ec  C Tst C 1 s es   es ¼ es  C st ec =C t  1  T dt ¼ dt þ eT a C s ea þ et et =C t  1  T det ¼ det þ eT C e þ e e =C a s t s c t C t ¼ C t  C Tst C 1 s C st

398

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402

the finite element method to seek a numerical solution. For the convenience of derivation, we introduce the following matrix notations:

 ¼ be11

V T Qw ¼ 0 with

2

T

2e12 e22 K 11 K 12 þ K 21 K 22 E2D1 E2D2 c

^ ¼ bw1 w2 w3 gcT w

ð41Þ

ð50Þ

Q ¼ hST Si

3 1 0 0 0 1 0 0 0 ... 60 1 0 0 0 1 0 0 ...7 6 7 wT ¼ k6 7 40 0 1 0 0 0 1 0 ...5 0

C ¼ bC11 2C12 C22 2C13 2C23 C33 E1 E2 E3 cT From Eqs. (11) and (14), we have the following

^ þ C  þ Cl 1 w ^ ;1 þ Cl2 w ^ ;2 C ¼ Ch w

ð42Þ

0 60 6 6 60 6 6 o 6 ox3 6 Ch ¼ 6 60 6 60 6 6 60 6 60 4 0 2 1 60 6 6 60 6 6 60 6 C ¼ 6 60 6 60 6 60 6 6 40 0

0

0

0

0

0

0 0

0 0

o ox3

0

0

o ox3

0

0

0

0

0

0

0 0 x3 1 0 0 1

0 0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

3

2

1 0 60 1 0 7 7 6 7 6 0 7 60 0 7 6 7 6 0 7 60 0 7 6 7 0 7 Cl1 ¼ 6 0 0 6 7 6 0 7 60 0 7 6 7 60 0 7 0 7 6 6 7 40 0 0 5 0 0  oxo3 3 0 0 0 0 x3 0 0 0 7 7 7 0 x3 0 0 7 7 7 0 0 0 07 7 0 0 0 07 7 Cl2 7 0 0 0 07 7 0 0 1 07 7 7 0 0 0 15 0 0 0 0

0

0

3

0 7 7 7 0 0 7 7 7 1 0 7 7 0 0 7 7 7 0 0 7 7 0 1 7 7 7 0 0 5 0 0 2 0 0 61 0 6 6 60 1 6 6 60 0 6 ¼6 60 0 6 60 0 6 60 0 6 6 40 0

1 2

Z *

"

CT

X

CE T

e

e C

d

# +

C dX

1 2

Z

0

ð51Þ

ð52Þ

The solution of warping can be expressed symbolically as

ð43Þ

b 0 V¼V

ð53Þ

Substituting Eq. (53) back into Eq. (48), one can obtain the total electric enthalpy density asymptotically correct up to the zeroth order of lg2 as

0 0 0 0 1 0 0 0 0

0

  b T Dh þ D  2H0 ¼ T V 0

3

7 7 7 7 7 7 0 7 7 0 7 7ð44Þ 7 0 7 7 0 7 7 7 1 5 0 0 0

hCT DCidX

ð45Þ

X

C

where C E , e, and d are the elastic material matrix, piezoelectric coefficient matrix, and dielectric coefficient matrix, respectively, with the partitioned matrices given in Eq. (31). The zeroth-order approximation of the electric enthalpy density in Eq. (34) can be written in the following matrix form:

H0 ¼

0 1 ...

FV ¼ ðQ wwT  IÞDh 

The electric enthalpy can be easily calculated by substituting Eqs. (42) into Eq. (2) as

H¼

0

where K is the Lagrange multiplier and can be obtained as K ¼ wT Dh  by taking advantage of the properties of the kernel matrix w. Substituting K back into Eq. (51), we obtain

0

0

0 1 0

FV þ Dh  ¼ Q wK

with the following operator matrices introduced:

2

0

where k can be determined if we normalize w such that wT Q w ¼ I with I as the identity matrix. Now our problem is transformed to find the stationary value of Eq. (48) numerically, subject to the constraints in Eq. (50). The Euler–Lagrange equation can be derived as

1 T h C D C0 i 2 0

ð46Þ

^ þ C . with C0 ¼ Ch w In order to solve the zeroth-order warping, we can discretize the transverse normal line into 1D finite elements and express the warping field as

^ i Þ ¼ Sðx3 ÞVðx1 ; x2 Þ wðx

ð47Þ

where S is the shape function and V is the nodal value of warping functions along the transverse normal. Substituting Eq. (47) into Eq. (46), one can obtain a discretized version of H0 as

2H0 ¼ V T FV þ 2V T Dh  þ T D 

ð48Þ

with

D E F ¼ h½Ch ST D½Ch Si Dh ¼ h½Ch ST DC i D ¼ CT DC

ð54Þ

which is an numerical version of Eq. (38). For construction of a refined model with better prediction, we need to carry out the asymptotical expansion of the electric enthalpy of the original 3D problem to a higher order. Based on the order assessment in Eq. (29) and VAM, the next approximation of the warping functions corresponds to the stationary point of the following functional

1 2



P ¼ hCT DCi  hPi wi i þ si wþi þ bi wi  Dþ dg þ  D dg 

ð55Þ

subjected to constraints in Eqs. (6) and (13). A discretized version of P can be derived as

2P ¼ V T FV þ 2V T ðDh  þ Dhl1 V ;1 þ Dhl2 V ;2 Þ þ T D  þ V T;1 Dl1 l1 V ;1 þ V T;2 Dl2 l2 V ;2   þ 2 V T;1 Dl1   þ V T;2 Dl2   þ V T;1 Dl1 l2 V ;2 þ 2V T L

ð56Þ þT

where L contains the load related terms such that L ¼ S s ST b  hST Pi with s ¼ bs1 s2 s3  Dþ cT , b ¼ bb1 b2 b3  D cT , P ¼ bP1 P2 P3 0cT , and

Dhl1 ¼ h½Ch ST D½Cl1 Si Dhl2 ¼ h½Ch ST D½Cl2 Si D ¼ hCT DC i Dl1 l1 ¼ h½Cl1 ST D½Cl1 Si Dl1 l2 ¼ h½Cl1 ST D½Cl2 Si Dl2 l2 ¼ h½Cl2 ST D½Cl2 Si

ð57Þ

Dl1  ¼ h½Cl1 ST DC i Dl2  ¼ h½Cl2 ST DC i We can simply perturb the warping function such that

V ¼ V0 þ V1

ð58Þ

with V 1 asymptotically smaller than V 0 . Substituting Eq. (58) into Eq. (56), one can obtain the leading terms for the first-order approximation as

2P1 ¼ V T1 EV 1 þ 2V T1 D1 ;1 þ 2V T2 D2 ;2 þ 2V T1 L

ð59Þ

where

ð49Þ

The discretized form of constraint in Eqs. (6) and (13) can be written as

  b 0  Dl  D1 ¼ Dhl1  DThl1 V 1   T b D2 ¼ Dhl2  Dhl2 V 0  Dl2 

ð60Þ ð61Þ

399

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402

It is understood that the order of the loads in Eq. (29) is associated with warping functions of different orders, as shown in [43]. Similarly as in the zeroth-order approximation, we can solve the first-order warping field as

V 1 ¼ V 11 ;1 þ V 12 ;2 þ V 1L

ð62Þ

Then we can obtain an approximation of the functional P in Eq. (55) asymptotically correct up to the order of lðh=lÞ2 g2 , given by

2P1 ¼ T A þ T;1 B;1 þ 2T;1 C ;2 þ T;2 D;2 þ 2T A;1 ð63Þ

where T b b T Dl l V B¼V 0 1 1 0 þ V 11 D1

  T T 1 b b T Dl l V C¼V 0 1 2 0 þ 2 V 11 D2 þ D1 V 12

T b b T Dl l V D¼V 0 2 2 0 þ V 12 D2

b 0 þ DT V b b 0 þ DT V b b T Dhl V b T Dhl V A¼V B¼V 0 0 l1  0 l2  0 1 2   b T L  1 DT V 1L;1 þ V T L;1 þ DT V 1L;2 þ V T L;2 C¼V 1 11 2 12 0 2 ð64Þ It is noted that the quadratic term involving applied loads is neglected because it will not affect the 2D plate model. 4. Transforming into a generalized Reissner–Mindlin model In order to obtain a refined 2D model that is of practical use, we can transform the asymptotically correct energy functional in Eq. (63) into a generalized Reissner–Mindlin model, in which the transverse shear strains are introduced as two additional degrees of freedom. Following [43], we can derive the following identity between the two sets of generalized strain measures:

 ¼ R  Da c;a

ð65Þ

with c ¼ b2c13 2c23 cT denoting the transverse shear strains, and

R¼



e11

2e12

e22

K 11

K 12

þ

K 21

K 22

E1

T E2

  D2 ¼

0 0 0 1 0 0 0 0

0 0 0

0 1 0 0 0

0 0 0 0 1 0 0 0

T T ð66Þ

0 0 0 0 0 1 0 0

With the identity in Eq. (65), we can express the enthalpy functional in Eq. (63) in terms of the generalized strain measures of the generalized Reissner–Mindlin model as T

2P1 ¼ R AR þ T

RT;1 BR;1

þ T

2RT;1 CR;2

þ

RT;2 DR;2

T

ð71Þ

F ¼ A  PT G1 DT1 A M ¼ C þ ðX  AÞDa G1 bm1 m2 cT;a If the model in Eq. (70) could be equivalent to the Reissner–Mindlin model in Eq. (68), the following equalities should hold:

X ¼ ðA þ PT G1 PÞ F R ¼ M 

RT;1 BR;1

L ¼ 0

2RT;1 CR;2

RT;2 DR;2

T

þ 2R AR;1 T

þ 2R BR;2 þ 2R C  2R AD1 c;1  2R AD2 c;2

L ¼ RT;1 BR;1 þ 2RT;1 CR;2 þ RT;2 DR;2

2PR ¼ RT XR þ 2cT PR þ cT Gc þ 2RT F R

ð68Þ

In order to find an equivalent Reissner–Mindlin model, we need to eliminate all the partial derivatives of the quantities in Eq. (67). We can derive the relation between shear strains and the derivative of R from the equilibrium equations balancing bending moments as

Gc þ PR ¼ DTa ðXR;a þ PT c;a þ F R;a Þ þ bm1 m2 cT

ð69Þ

By taking advantage of Eq. (69), Eq. (67) can be rewritten as

2P1 ¼ RT ðA þ PT G1 PÞR þ cT Gc þ 2cT PR þ 2RT FR;1 þ 2RT JR;2 þ RT;1 BR;1 þ 2RT;1 CR;2 þ RT;2 DR;2 þ 2RT M

ð70Þ

T

ð73Þ

with

B ¼ B þ AD1 G1 DT1 A 1

D ¼ D þ AD2 G

ð74Þ

DT2 A

Now, driving L to be zero is greatly simplified, which can be achieved by some optimization techniques. Having constructed the generalized Reissner–Mindlin model, we can then approximate the original electromechanical problem governed by the 3D Hamilton’s extended principle in Eq. (27) using the following 2D variational statement:

Z

t2

½dðK2D  Pg Þ þ dW2D dt ¼ 0

ð75Þ

t1

with K2D obtained from Eq. (17), dW2D from Eq. (24), and Pg as the 2D electric enthalpy defined for the generalized Reissner–Mindlin model as

ð67Þ

For a fully-coupled electromechanical problem, a generalized Reissner–Mindlin model can be of the following form

ð72Þ T

with L ¼ þ þ þ 2R FR;1 þ 2R JR;2 . Here X; P, G, and F R are what we need for the construction of the generalized Reissner–Mindlin model. For general cases, it is not possible for L to be zero and the best one can do is to use some optimization techniques to derive L to be as close to zero as possible. If we restrict the material to be monoclinic with in-plane polarization, which is common for practical applications of smart plates for in-plane sensing and actuation, the stiffness matrices A and B vanish in Eq. (67). Therefore, we don’t have first-order terms in Eq. (67). Consequently, the matrix related with the first-order terms in the generalized Reissner–Mindlin Model, P, is zero, X is identically equal to A, and L is simplified to be the following quadratic form

C ¼ C þ AD1 G1 DT2 A

denoting the generalized strains defined in the generalized Reissner–Mindlin model which are the counterparts of  defined previously, and

D1 ¼

B ¼ B þ XD1 G1 DT1 ð2A  XÞ i 1h C ¼ C þ ð2A  XÞD1 G1 DT2 X þ XD1 G1 DT2 ð2A  XÞ 2 D ¼ D þ XD2 G1 DT2 ð2A  XÞ J ¼ B  PT G1 DT2 A

þ 2T B;2 þ 2T C b T Dh þ D A¼V 0

where

Pg ¼

Z

Pf dX

ð76Þ

X

where

2Pf ¼ RT AR þ cT Gc þ 2RT F R

ð77Þ

5. Numerical examples The finite element formulation has been implemented into the variational asymptotic plate and shell analysis (VAPAS), a generalpurpose 1D through-the-thickness analysis code for plate and shell modeling. With the electric potential introduced as an independent degree of freedom, VAPAS is now able to model composite plates with embedded or attached piezoelectric material. To demonstrate the accuracy and efficiency of the present theory, a few examples will be analyzed using VAPAS and the results will be

400

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402

compared with those available in the literature or the original 3D model. First, we evaluate the stiffnesses needed for the coupled 2D model. A piezoelectric composite plate in cylindrical bending is used as an example. The length of the plate along x2 is 0.1 m, and it is infinite along x1 . The thickness of the plate is 0.01 m. This plate is made of PZT-4 with the material properties provided in Table 2. It should be noted that the material properties are taken from Heyliger’s paper [2] (see Table 2) in which the properties are used for the polarization along the thickness direction of the plate ðx3 Þ. If the piezoelectric material is poled along an in-plane direction ðx2 Þ, the constitutive relations of stresses, strains, electric displacements, and electric fields can be obtained from a 90° rotation around the other in-plane direction ðx1 Þ followed by a 180° rotation around the thickness direction ðx3 Þ. For this single layer plate, we consider the cases of in-plane polarization and thickness polarization. For the case of in-plane polarization, the poling direction of the piezoelectric materials is assumed to be in x2 direction. These stiffnesses are calculated using VAPAS and all the nonzero values are listed compared with Ref. [44] in Tables 3 and 4. It should be noted that Ref. [44] is only applicable to single layer plate, and our model can deal with the multi-layer problems. It is found out that for the case of in-plane polarization (see Table 3), the stiffness terms A11 , A13 , A22 , A33 , A44 , A46 , A55 , and A66 are the same as those obtained using the classical lamination plate theory (CLPT). In addition, the present theory predicts stiffness related to electric fields A18 , A27 , A38 , A77 , and A88 . It can be seen that Af 7 and Af 8 with f ranging from four to six are all zero, which means there is no coupling between in-plane electric fields and curvatures. As shown in Table 4, for the case of thickness polarization, the stiffnesses A11 , A13 , A33 , A44 , A46 , and A66 are larger than those of CLPT and the differences depend on piezoelectric and dielectric properties of the material. There is no electromechanical coupling in this material polarization case. Next, we calculate the displacements of this single layer plate with prescribed electric potential on the lateral boundary. We con-

Table 2 Material properties of piezoceramic materials

Table 4 Comparison of stiffnesses for thickness polarization Stiffness 9

A11 ð10 N=mÞ A13 ð108 N=mÞ A22 ð108 N=mÞ A33 ð109 N=mÞ A44 ð103 NmÞ A46 ð103 NmÞ A55 ð103 NmÞ A66 ð103 NmÞ A77 ð1010 C=VÞ A88 ð1010 C=VÞ

Le

Present

CLPT

1.07663 4.64896 3.06000 1.07663 8.97195 3.87414 2.55000 8.97195 1.93740 1.93740

1.07663 4.64896 3.06000 1.07663 8.97195 3.87414 2.55000 8.97195 1.93740 1.93740

0.911681 2.99943 3.06000 0.911681 7.59734 2.49953 2.55000 7.59734 N/A N/A

sider the applied electric potential at x2 ¼ 0 and x2 ¼ 0:1 m are 0 V and 1.0 V, respectively. For this plate, the 2D solution can be obtained from Ritz approximation method. Let us assume that the 2D field variables are uniform along the x1 direction, and their distribution along x2 can be assumed to be the following forms:

u1 ðx2 Þ ¼ U 1 cos



p

Lt 

 x2

 x2 Lt   p x2 u3 ðx2 Þ ¼ U 3 sin Lt   p u1 ðx2 Þ ¼ u1 cos x2 Lt   p u2 ðx2 Þ ¼ u2 cos x2 Lt u2 ðx2 Þ ¼

U 2

cos

p

ð78Þ

Uðx2 Þ ¼ a1 x22 þ a2 x2 þ a3 in which u1 , u2 , and u3 are 2D displacements, u1 and u2 are rotations around x2 and x1 , U is the 2D potential, Lt is length of the plate, and the three unknowns a1 , a2 , and a3 can be reduced to one by applying two boundary conditions of prescribed potential. Substituting the expressions of displacements and electric potential in Eq. (78) into the expression of generalized strains of Reissner–Mindlin model, we have T

Quantity

PZT-4

Quantity

PZT-4

E1 ðGPaÞ E2 ðGPaÞ E3 ðGPaÞ G23 ðGPaÞ G13 ðGPaÞ G12 ðGPaÞ d11 =0 d22 =0

81.3 81.3 64.5 25.6 25.6 30.6 1475 1475

m12 m13 m23

0.329 0.432 0.432 12.72 5.2 5.2 15.08 1300

e24 ðC=m2 Þ e31 ðC=m2 Þ e32 ðC=m2 Þ e33 ðC=m2 Þ d33 =0

Table 3 Comparison of stiffnesses for in-plane polarization Stiffness

Le

Present

CLPT

A11 ð109 N=mÞ A13 ð108 N=mÞ A22 ð108 N=mÞ A33 ð109 N=mÞ A44 ð103 NmÞ A46 ð103 NmÞ A55 ð103 NmÞ A66 ð103 NmÞ A18 ð102 C=mÞ A27 ðC=mÞ A38 ðC=mÞ A77 ð1010 C=VÞ A88 ð1010 C=VÞ G11 ð108 N=mÞ G22 ð108 N=mÞ

9.54292 3.27065 2.56000 7.57095 7.95243 2.72554 2.13333 6.30913 2.05278 0.12720 0.18405 1.30538 1.16556 N/A N/A

9.54292 3.27065 2.56000 7.57095 7.95243 2.72554 2.13333 6.30913 2.05278 0.12720 0.18405 1.30538 1.16556 2.76730 2.31848

9.54292 3.27065 2.56000 7.57095 7.95243 2.72554 2.13333 6.30913 N/A N/A N/A N/A N/A N/A N/A

R ¼ ½ e11 2e12 e22 K 11 K 12 þ K 21 K 22 E1 E2  h iT ou 1 ou 1 ou1 ou2 oU oU 1 2 þ ou þ ooxu12 ooxu22  ox  ¼ ou ð79Þ ox1 ox2 ox1 ox2 ox1 ox2 ox 1 2 h iT 3 3 u2 þ ou c ¼ u1 þ ou ð80Þ ox1 ox2 For this problem, the generalized force in Reissner–Mindlin model F R is zero. Substituting Eqs. (79) and (80) into Eq. (77), we can obtain the energy functional Pf . The six unknowns U 1 , U 2 , U 3 , u1 , u2 , and one coefficient in U can be solved by applying Ritz method and taking the derivatives of Pf with respect to these six unknown coefficients. The electric potential turns out to be a linear function of x2 , which accurately predicts the real distribution of original 3D structures as

U ¼ /1 þ

x2 ð/2  /1 Þ Lt

ð81Þ

where /1 and /2 are the prescribed potential at x2 ¼ 0, x2 ¼ Lt , respectively. The displacements of our model are compared with Table 5 Displacements of the single layer PZT-4 plate Displacements

3D

2D

Difference (%)

U 2mid ð1011 mÞ U 3min ð1011 mÞ U 3max ð1011 mÞ

9.79150 1.17391 1.17391

9.85251 1.13437 1.13437

0.62 3.37 3.37

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402 Table 6 Stiffnesses for the two layer plate Stiffness 9

A11 ð10 N=mÞ A13 ð108 N=mÞ A22 ð108 N=mÞ A33 ð109 N=mÞ A44 ð103 NmÞ A46 ð103 NmÞ A55 ð103 NmÞ A66 ð103 NmÞ A18 ð102 C=mÞ A27 ðC=mÞ A38 ðC=mÞ A77 ð1010 C=VÞ A88 ð1010 C=VÞ G11 ð108 N=mÞ G12 ð106 N=mÞ G22 ð108 N=mÞ

9.54292 3.27065 2.56000 7.57095 7.95243 2.72554 2.13333 6.30913 2.05278 0.12720 0.18405 1.30538 1.16556 2.17017 7.96231 1.87486

Table 7 Displacements of the two layer PZT-4 plate Displacements

3D

2D

Difference (%)

U 1min ð1011 mÞ U 1max ð1011 mÞ U 2min ð1011 mÞ U 2max ð1010 mÞ U 3min ð1011 mÞ U 3max ð1011 mÞ

4.23167 4.14408 8.51073 1.00459 3.75704 5.97412

4.18757 4.18757 8.56344 1.01077 3.81259 5.95582

1.04 1.05 0.62 0.62 1.48 0.31

those of 3D solution, which is presented in Table 5. The 3D solution is also obtained from Ritz Method. The next example is a two-layer piezoelectric composite plate under cylindrical bending. The plate is made of PZT-4 with a layup of [0/30]. The length of the plate is 0.1 m and the thickness is 0.01 m. The stiffness of this plate is given in Table 6. The comparison of minimum and maximum displacements are provided in Table 7. It is shown that our 2D models achieve an excellent accuracy in comparison with 3D solution. It should be noted that the Ritz method is approximate and its accuracy depends on the approximation and given conditions. For multi-layer plates and complex conditions of geometry, loading, and boundary, we need to resort to numerical methods to obtain 2D solution. One can use the 2D energy in Eq. (75) to develop a 2D coupled solver in terms of R and c to calculate the global responses of laminated piezoelectric plates similarly as what has been done for mechanical field in Ref. [45]. The 2D plate analysis based on the variational statement of Eq. (75) has 10 generalized strains (eight in R and two transverse shear strains), which can be derived from five mechanical degrees of freedom (three displacements of the reference plane plus two rotations of the transverse normal) and one electric degree of freedom (the 2D electric potential). For general plate analysis, we need to rely on plate finite elements to find numerical solution for the 2D variational statement in Eq. (75), which means we need to augment the regular Reissner–Mindlin plate element with one more nodal degree of freedom to account for the 2D electric potential. How to develop the generalized Reissner–Mindlin plate elements capable of dealing with laminated piezoelectric plates will be studied in future work. 6. Conclusion In this paper, a unique electromechanical Reissner–Mindlin model has been constructed for piezoelectric composite plates without any prescribed electric potential through the thickness. Four fundamental variables including three displacements and

401

electric potential, which are independent of laminate layer number, are used in developing the proposed theory. VAM is applied to perform the dimensional reduction from 3D problem to a 1D through-the-thickness analysis and a 2D plate analysis. Shear deformation effect has been introduced to construct the refined model and shear stiffness is solved by minimization techniques. With the solved warping functions, the total energy asymptotically correct up to the second order can be expressed in terms of 2D strains, curvature, and electric field variables. The accuracy and application of this model have been demonstrated using a few examples. Acknowledgement The present work is supported, in part, by the National Science Foundation under Grant DMI-0522908 and the Space Dynamics Laboratory under Shunk Works Grant and Enabling Technologies Grant. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsement of the funding agencies. References [1] Ray MC, Rao KM, Samanta B. Exact analysis of coupled electroelastic behavior of a piezoelectric plate under cylindrical bending. Comput Struct 1992;45(4):667–77. [2] Heyliger PR, Brooks S. Exact solutions for laminated piezoelectric plates in cylindrical bending. J Appl Mech 1996;63:903–10. [3] Vel SS, Batra RC. Exact solution for rectangular sandwich plates with embedded piezoelectric shear actuators. AIAA J 2001;39(7):1363–73. [4] Allik H, Hughes TJR. Finite element method for piezoelectric vibration. Int J Numer Meth Eng 1970;2:151–7. [5] Dimitriadis EK, Fuller CR, Rogers CA. Piezoelectric actuators for distributed vibration excitation of thin plates. J Vib Acoust 1991;113:100–7. [6] Crawley EF, Lazarus KB. Induced strain actuation of isotropic and anisotropic plates. AIAA J 1991;29(6):944–51. [7] Wang BT, Rogers CA. Laminate plate theory for spatially distributed induced strain actuators. J Compos Mater 1991;25:433–52. [8] Hong CH, Chopra I. Modeling and validation of induced strain actuation of composite coupled plates. AIAA J 1999;37(3):372–7. [9] Ha SK, Keilers C, Chang FK. Finite element analysis of composite structures containing distributed piezoceramic sensors and actuators. AIAA J 1992;30(3):772–80. [10] Ray MC, Bhattacharya R, Samanta B. Static analysis of an intelligent structure by the finite-element method. Comput Struct 1994;52(4):617–31. [11] Detwiler DT, Shen MHH, Venkayya VB. Finite-element analysis of laminated composite structures containing distributed piezoelectric actuators and sensors. Finite Element Anal Design 1995;20(2):87–100. [12] Suleman A, Venkayya VB. A simple finite element formulation for a laminated composite plate with piezoelectric layers. J Intell Mater Syst Struct 1995;6(6):776–82. [13] Yu JY, Kang WY, Kim SJ. Elastic tailoring of laminated composite plate by anisotropic piezoelectric polymers – theory, computation, and experiment. J Compos Mater 1995;29(9):1201–21. [14] Wang J, Yong YK, Imai T. Finite element analysis of the piezoelectric vibrations of quartz plate resonators with higher-order plate theory. Int J Solids Struct 1999;36(15):2303–19. [15] Krommer M, Irschik H. A Reissner–Mindlin-type plate theory including the direct piezoelectric and the pyroelectric effect. Acta Mech 2000;141:51–69. [16] Wang J, Yang J. Higher-order theories of piezoelectric plates and applications. Appl Mech Rev 2000;53(4):87–99. [17] Kapuria S, Dube GP, Dumir PC. First-order shear deformation theory solution for a circular piezoelectric composite plate under axisymmetric load. Smart Mater Struct 2003;12:417–23. [18] Pai PF, Nayfeh AH, Oh K, Mook DT. Refined nonlinear model of composite plates with integrated piezoelectric actuators and sensors. Int J Solids Struct 1993;30(12):1603–30. [19] Icardi U, Sciuva MD. Large-deflection and stress analysis of multilayered plates with induced-strain actuators. Smart Mater Struct 1995;5:140–64. [20] Saravanos D, Heyliger P, Hopkins D. Layerwise mechanics and finite element for the dynamic analysis of piezoelectric composite plates. Int J Solids Struct 1997;34(3):359–78. [21] Zhou X, Chattopadhyay A, Thornburgh R. Analysis of piezoelectric smart composites using a coupled piezoelectric-mechanical model. J Intell Mater Syst Struct 2000;11:169–79. [22] Sheikh AH, Topdar P, Halder S. An appropriate FE model for through-thickness variation of displacement and potential in thin/moderately thick smart laminates. Compos Struct 2001;51(4):401–9.

402

L. Liao, W. Yu / Composite Structures 88 (2009) 394–402

[23] Cen S, Soh A, et al. A new 4-node quadrilateral FE model with variable electric degrees of freedom for the analysis of piezoelectric laminated composite plates. Compos Struct 2002;58(4):583–99. [24] Kapuria S, Achary GGS. A coupled consistent third-order theory for hybrid piezoelectric plates. Compos Struct 2005;70(1):120–33. [25] Kapuria S, Achary GGS. Nonlinear coupled zigzag theory for buckling of hybrid piezoelectric plates. Compos Struct 2006;74(3):253–64. [26] Zhen W, Wanji C. Refined triangular element for laminated elasticpiezoelectric plates. Compos Struct 2007;78(1):129–39. [27] Chee C, Tong L, Steven GP. Piezoelectric actuator orientation optimization for static shape control of composite plates. Compos Struct 2002;55(2):169–84. [28] Nguyen Q, Tong L. Shape control of smart composite plate with nonrectangular piezoelectric actuators. Compos Struct 2004;66(1–4):207–14. [29] Zhang HY, Shen YP. Vibration suppression of laminated plates with 1–3 piezoelectric fiber-reinforced composite layers equipped with interdigitated electrodes. Compos Struct 2007;79(2):220–8. [30] Liao L, Yu W. Asymptotical construction of a fully coupled, Reissner–Mindlin model for piezoelectric composite plates. Smart Mater Struct 2007;17(1):1–14. [31] Warner AWJ, Goldfrank B. Lateral field resonators. In: Proceedings of the annual frequency control symposium; 1985. p. 473–4. [32] Abe H, Yoshida T, Watanabe H. Energy trapping of thickness-shear vibrations excited by parallel electric field and its application to piezoelectric vibratory gyroscopes. In: Proceedings of the IEEE ultrasonics symposium, vol. 1, 1998. p. 467–72. [33] Khan A, Ballato A. Piezoelectric coupling factor calculations for plates of langatate driven in simple thickness modes by lateral-field-excitation. IEEE Trans Ultrasonics Ferroelectrics Frequency Control 2002;49(7):922–8.

[34] Nedorezov S, Shmaliy O, Shmaliy Y. On influence of an electrode edge potential upon frequencies of life piezoelectric resonators employing thickness-shear vibrations. In: Proceedings of the 2004 IEEE international frequency control symposium and exposition, vol. 1, 2004. p. 585–90. [35] Yamada T, Niizeki N. Formulation of admittance for parallel field excitation of piezoelectric plates. J Appl Phys 1970;41(9):3604–9. [36] Ikeda T. Depolarizing-field effect in piezoelectric resonators. Ferroelectrics 1982;43(1–2):3–15. [37] Zhao Z. Universal resonant-frequency equation for the thickness vibrations of piezoelectric plate. Shengxue Xuebao/Acta Acustica 1996;21(4):332–8. [38] Zhao Z. Further study on the thickness vibration of piezoelectric plate excited by parallel field. Shengxue Xuebao/Acta Acustica 1996;21(2):156–64. [39] Lee PCY. Electromagnetic radiation from an at-cut quartz plate under lateralfield excitation. In: Ultrasonics Symposium Proceedings; 1988. p. 407–11. [40] Berdichevsky VL. Variational-asymptotic method of constructing a theory of shells. PMM 1979;43(4):664–87. [41] Yang J. The mechanics of piezoelectric structures. 1st ed. World Scientific Publishing Company; 2006. [42] Danielson DA. Finite rotation with small strain in beams and plates. In: Proceedings of the 2nd Pan American congress of applied mechanics, Valparaiso, Chile; 1991. [43] Yu W, Hodges DH, Volovoi VV. Asymptotic construction of Reissner-like models for composite plates with accurate strain recovery. Int J Solids Struct 2002;39(20):5185–203. [44] Le KC. Vibrations of shells and rods. 1st ed. Germany: Springer; 1999. [45] Hodges DH, Atılgan AR, Danielson DA. A geometrically nonlinear theory of elastic plates. J Appl Mech 1993;60(1):109–16.