The exact solution of coupled thermoelectroelastic behavior of piezoelectric laminates

Computers and Structures 80 (2002) 1201–1212 www.elsevier.com/locate/compstruc

The exact solution of coupled thermoelectroelastic behavior of piezoelectric laminates C. Zhang

a,* ,

Y.K. Cheung b, S. Di a, N. Zhang

a

a

b

Faculty of Engineering, University of Technology, Sydney, Broadway 1, NSW 2007, Australia Department of Civil and Structural Engineering, University of Hong Kong, Pokfulam Road, Hong Kong Received 24 July 2001; accepted 15 February 2002

Abstract Exact solutions for static analysis of thermoelectroelastic laminated plates are presented. In this analysis, a new concise procedure for the analytical solution of composite laminated plates with piezoelectric layers is developed. A simple eigenvalue formula in real number form is directly developed from the basic coupled piezoelectric differential equations and the difficulty of treating imaginary eigenvalues is avoided. The solution is defined in the trigonometric series and can be applied to thin and thick plates. Numerical studies are conducted on a five-layer piezoelectric plate and the complexity of stresses and deformations under combined loading is illustrated. The results presented here could be used as a benchmark for assessing any numerical solution by approximate approaches such as the finite element method while also providing useful physical insight into the behavior of piezoelectric plates in thermal environment. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Thermoelectroelastic; Laminate; Composite; Piezoelectric; Exact solution; Trigonometric series

1. Introduction There has been a great deal of interest recently in the study of smart structures. With its monitoring and adaptive controlling capacities, a smart structure system has found its wide applications in space structures, aerocraft, automobile and medical equipment. The piezoelectric materials are used in a smart structure system because of its properties for inducing direct and converse piezoelectric effects. The coupled effects of mechanic and electric properties of piezoelectric materials make structural analyses complicated. Piezoelectric plates composed of orthotropic or anisotropic layers are usually analyzed by using the laminated composite plate theory. For the analytical solution of laminated composite plates, Pagano [1] firstly presented the exact solution of the laminated plate under cylindrical bending. Ray et al. [2–4], Heyliger and Brooks [5,6] extended this methodology to develop the exact solution of laminated plates

*

Corresponding author.

composed either entirely or in part of piezoelectric layers. Senthil and Batra [7] recently presented the analytical solutions of piezoelectric laminated plates via Eshelby–Stroh formalism. Further, taking the thermal effect into account, Xu et al. [8] presented a threedimensional analysis for the coupled thermoelectroelastic response of mutilayered hybrid composite plates with four simply supported edges by using the state space method. Dube [9] developed a series solution for thermoelectroelastic plate with single layer under cylindrical bending by using Pagano’s approach. Meanwhile, a number of papers, see Sung and Charles [10], Hwang and Park [11], and Chattopadhyay et al. [12], have reported the effect of thermal loading to smart piezoelectric structures by using the finite element method. In the above mentioned analytical solutions to laminated piezoelectric plates, complex mathematical derivations (see, for instance Refs. [5,8]) are needed in order to circumvent the difficulty of treating imaginary eigenvalues and obtaining real eigensolutions of the differential equations. In this paper, a concise solution procedure is directly developed from the basic differential equations

0045-7949/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 2 ) 0 0 0 6 0 - 3

1202

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

and used to obtain the analytical solutions of thermopizoelectroelastic laminated plates. The new solution procedure could be considered as an alternative way of exact solutions to piezoelectric laminated plates. Moreover, with a number of piezoelectric finite element models [13] to cope with complicated boundary problems of smart structures, it is very important to develop a simple procedure to obtain analytical solutions in a few simple geometries for comparison. In this analysis, considering the simply supported, grounded and laminated rectangular plates with piezoelectric layers, field variables in basic equations are expanded in the trigonometric series and the unknown function amplitudes are determined by the continuous interface and the surface conditions. The numerical example includes different loading cases and the presented results provide a better understanding of the interaction between the mechanical, thermal and electric fields and also may be used as a benchmark for some approximate piezoelectric plate theories, finite element and boundary element methods.

2. Mathematical formulation Fig. 1 shows the geometrical configuration of a piezoelectric laminated plate with simply supported and grounded edges. The thickness dimension, H, of the laminated plate coincides with the z-direction. The longitudinal length, a, and the width, b, of the plate are along the x-direction and y-direction respectively. The plate consists of arbitrary n layers with the top and bottom of piezoelectric material serving as an actuator or a sensor. The analytical formulation is based on the

linear three-dimensional anisotropic theory with the consideration of thermal effects. The linear constitutive equations for any layer (z is considered as local coordinate,
hk 6 z 6 hk ) are ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðk Þ ðkÞ ðkÞ ðkÞ ðkÞ rðkÞ ; x ¼ C11 u;x þ C12 v;y þ C13 w;z
e31 Ez
k1 T

ð1Þ ðkÞ

ðkÞ

ðk Þ

ðkÞ

ðkÞ

ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ rðkÞ ; y ¼ C12 u;x þ C22 v;y þ C23 w;z
e32 Ez
k2 T

ð2Þ ðkÞ

ðkÞ

ðk Þ

ðkÞ

ðkÞ

ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ; rðkÞ z ¼ C13 u;x þ C23 v;y þ C33 w;z
e33 Ez
k3 T

ð3Þ
 ðkÞ ðkÞ ðkÞ
e24 Eyk ; sðkÞ yz ¼ C44 v;z þ w;y

ð4Þ


 ðkÞ ðkÞ þ v ; ¼ C u sðkÞ 66 xy ;y ;x

ð5Þ


 ðkÞ ðkÞ ðkÞ
e15 Exk ; sðkÞ xz ¼ C55 u;z þ w;x

ð6Þ


 ðkÞ ðkÞ ðk Þ ðkÞ þ j11 ExðkÞ ; DðkÞ x ¼ e15 u;z þ w;x

ð7Þ


 ðkÞ ðk Þ ðk Þ ðkÞ þ j22 EyðkÞ ; DðkÞ y ¼ e24 v;z þ w;y

ð8Þ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðk Þ ðkÞ ðkÞ ðkÞ ðkÞ DðkÞ ; z ¼ e31 u;x þ e32 v;y þ e33 w;z þ j33 Ez þ r3 T

ð9Þ where (k) denotes the order of layer and the suffix 1, 2, 3 refer to x-, y-, z-directions; u, v, w are the displacement components along the direction of coordinate x, y, z; Ex , Ex , Ez are electric field components; T is temperature change and Dx , Dy , Dz are electrical displacement; Cij and eij are the elastic and piezoelectric coefficient respectively; ri and ki are pyroelectric coefficients and stress–temperature coefficients; j11 , j22 , j33 are permittivities. The subscript comma represents the partial differentiation with respect to the following quantitative. The stress equilibrium equations are ðkÞ ðkÞ rðkÞ x;x þ sxy;y þ sxz;z ¼ 0;

ð10Þ

ðkÞ ðkÞ sðkÞ xy;x þ sy;y þ ryz;z ¼ 0;

ð11Þ

ðkÞ ðkÞ sðkÞ xz;x þ syz;y þ rx;z ¼ 0;

ð12Þ

where the inertial terms with respect to time are not included under consideration of static analysis. The charge equilibrium is Fig. 1. Geometry of the composite laminates with piezoelectric layers.

ðkÞ ðkÞ DðkÞ x;x þ Dy;y þ Dz;z ¼ 0:

ð13Þ

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212 ðkÞ

pzðkÞ ¼
k33 T;zðkÞ

The thermal conduction equilibrium is px;x þ py;y þ pz;z ¼ 0;

ð14Þ

where px , py , pz are heat flux components in the directions of x-, y- and z-axis. The relations between the electric field components and electrostatic potential are ExðkÞ ¼


o/ðkÞ ; ox

ð15Þ

EyðkÞ ¼


o/ðkÞ ; oy

ð16Þ

EzðkÞ ¼


o/ðkÞ : oz

ð17Þ

The relations between the temperature change and the heat flux components are ¼

ðkÞ
k11 T;xðkÞ ;

pyðkÞ ¼

ðkÞ
k22 T;yðkÞ ;

pxðkÞ

ð18Þ ð19Þ

1203

ð20Þ

where / is the electric potential and k11 , k22 , k33 are the thermal conductivity coefficients referring to three directions and all the permittivities and pyroelectric coefficients are zero in the case of the non-piezoelectric substrate. In the analysis of composite laminates, the displacement, electric potential and transverse stresses, electric displacement and temperature change which are required to satisfy the interface continuities through the thickness direction, form the basic field variable space for the governing equations. So the interface continuous conditions for the thermoelectroelastic plate are     hk hkþ1 uðkÞ x; y; ¼ uðkþ1Þ x; y;
; ð21Þ 2 2     hk hkþ1 ¼ vðkþ1Þ x; y;
; vðkÞ x; y; 2 2 ðkÞ

w



hk x; y; 2



ðkþ1Þ

¼w



 hkþ1 ; x; y;
2

ð22Þ

ð23Þ

Fig. 2. sxz through of the plate in loading: (a) case 1 (x ¼ 0, y ¼ b=2), (b) case 2 (x ¼ 0, y ¼ b=2) and (c) case 3 (x ¼ 0, y ¼ b=2).

1204

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

Fig. 3. syz through thickness of the plate in loading: (a) case 1 (x ¼ a=2, y ¼ 0), (b) case 2 (x ¼ a=2, y ¼ 0) and (c) case 2 (x ¼ a=2, y ¼ 0).

rðkÞ z

sðkÞ xz







hk x; y; 2

hk x; y; 2

¼

rðkþ1Þ z

¼

sðkþ1Þ yz



 hkþ1 ; x; y;
2

ð24Þ

 hkþ1 ; x; y;
2

ð25Þ





    hk hkþ1 ðkþ1Þ ¼ s ; x; y; x; y;
sðkÞ yx yz 2 2     hk hkþ1 ¼ T ðkþ1Þ x; y;
; T ðkÞ x; y; 2 2     hk hkþ1 ¼ pzðkþ1Þ x; y;
; pzðkÞ x; y; 2 2     hk hkþ1 ¼ /ðkþ1Þ x; y;
; /ðkÞ x; y; 2 2 DðkÞ z



hk x; y; 2

 ¼

Dðkþ1Þ z



 hkþ1 : x; y;
2

and bottom surfaces of the plate, are functions of x and y. These prescribed quantities are expanded in the form of the trigonometric series as follows:
qiz ð x; y Þ ¼

X


qizmn sin am x sin bn y;

ð31Þ

m;n

ð26Þ


qixz ð x; y Þ ¼

X


qixzmn cos am x sin bn y;

ð32Þ


qiyzmn sin am x cos bn y;

ð33Þ

/
imn sin am x sin bn y;

ð34Þ


imn sin am x sin bn y; D

ð35Þ

i T
mn sin am x sin bn y;

ð36Þ

m;n

ð27Þ


qiyz ð x; y Þ ¼

X m;n

ð28Þ

/
i ð x; y Þ ¼

X m;n

ð29Þ

or
i ð x; y Þ ¼ D

ð30Þ

It is assumed that the traction components, the temperature variation and electric potential at the top

X m;n

T
i ð x; y Þ ¼

X m;n

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212


i ðx; yÞ; Dz ð x; y;  H Þ ¼ D z

or
pzi ð x; y Þ ¼

X

i
pmn sin am x sin bn y;

ð37Þ

m;n

where i ¼ 1,2 correspond to the top and bottom surfaces and the upper line of the letters means that the designated variables are known. On the edges of the plate, the displacements and potential conditions are assumed as

or

wðkÞ ¼ vðkÞ ¼ /ðkÞ ¼ 0;

ðx ¼ 0; aÞ;

ð38Þ

pz ð x; y;  H Þ ¼
pzi ðx; yÞ:

wðkÞ ¼ uðkÞ ¼ /ðkÞ ¼ 0;

ðy ¼ 0; bÞ:

ð39Þ

3. Solution of equations

rz ð x; y;  H Þ ¼
qiz ðx; yÞ;

ð40Þ

sxz ð x; y;  H Þ ¼
qixz ðx; yÞ;

ð41Þ

ð43Þ

or /ð x; y;  H Þ ¼ /
i ð x; y Þ;

The boundary conditions on the top and bottom are mechanical conditions

1205

ð44Þ

thermal conditions T ð x; y;  H Þ ¼ T
i ðx; yÞ;

ð45Þ

ð46Þ

The displacement, potential and temperature functions as well as the out-of-plane stresses satisfying Eqs. (38) and (39) may be expanded in double trigonometric series as follows: X ðkÞ f1mn ðzÞ cos am x sin bn y; ð47Þ uðkÞ ¼ m;n

syz ð x; y;  H Þ ¼


qiyz ðx; yÞ;

electric conditions

ð42Þ vðkÞ ¼

X

ðkÞ

f2mn ðzÞ sin am x cos bn y;

ð48Þ

m;n

Fig. 4. u through thickness of the plate in loading: (a) case 1 (x ¼ 0, y ¼ b=2), (b) case 2 (x ¼ 0, y ¼ b=2) and (c) case 3 (x ¼ 0, y ¼ b=2).

1206

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

Fig. 5. v through thickness of the plate in loading: (a) case 1 (x ¼ a=2, y ¼ 0), (b) case 2 (x ¼ a=2, y ¼ 0) and (c) case 3 (x ¼ a=2, y ¼ 0).

X

wðkÞ ¼

ðkÞ

f3mn ðzÞ sin am x sin bn y;

ð49Þ

pzðkÞ ¼

m;n

ðkÞ

f4mn ðzÞ sin am x sin bn y;

ð50Þ

m;n

sðkÞ xz ¼

X

ðkÞ

ð51Þ

ðkÞ

ð52Þ

f5mn ðzÞ cos am x sin bn y;

m;n

sðkÞ yz ¼

X

f6mn ðzÞ sin am x cos bn y;

m;n

rðkÞ z ¼

X

ðkÞ

f7mn ðzÞ sin am x sin bn y;

ð53Þ

m;n

DðkÞ z ¼

ðkÞ

f10mn ðzÞ sin am x sin bn y:

ð56Þ

mn

X

/ðkÞ ¼

X

X

ðkÞ

f8mn ðzÞ sin am x sin bn y;

ð54Þ

The continuity of displacements and out-of-plane stresses are realized by the amplitude functions about z-coordinate of Eqs. (47)–(56) satisfying the interface continuous conditions Eqs. (21)–(30). From Eqs. 14, ðkÞ ðkÞ (18)–(20), f9mn ðzÞ and f10mn ðzÞ can be solved independently. Considering simultaneously Eqs. (14), (18)–(20), (27), (28), (45), (46) and Eqs. (55) and (56), we obtain ðkÞ

ðkÞ

f9mn;;z ðzÞ
cðkÞ mn f9mn ðzÞ ¼ 0; h
. i1=2 ðkÞ ðkÞ ðkÞ k33 where, cðkÞ a2m k11 þ b2n k22 . mn ¼

ð57Þ

Eq. (57) is a set of second-order homogeneous linear differential equations and the solution is in form of ðkÞ

ðkÞ

ðkÞ


cmn z cmn z þ BðkÞ ; f9mn ðzÞ ¼ AðkÞ mn e mn e

ð58Þ

m;n ðkÞ

T

ðkÞ

¼

X m;n

ðkÞ f9mn ðzÞ sin am x sin bn y;

ð55Þ

ðkÞ ðkÞ

f10mn ðzÞ ¼
k33 f9mn;z ðzÞ
 ðkÞ ðkÞ ðkÞ ðkÞ
cmn z cmn z
BðkÞ ; ¼ cðkÞ mn k33 Amn e mn e

ð59Þ

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

1207

Fig. 6. w through thickness of the plate in loading: (a) case 1 (x ¼ a=2, y ¼ b=2), (b) case 2 (x ¼ a=2, y ¼ b=2) and (c) case 3 (x ¼ a=2, y ¼ b=2).

ðkÞ where the unknowns, AðkÞ mn , Bmn , are determined by the continuous conditions on the layer interfaces Eqs. (27) and (28), and the prescribed top and bottom surface conditions Eqs. (45) and (46). Using the solved T ðkÞ and pzðkÞ and setting the in-plane stresses rx , ry , and sxy , the in-plane electric displacements Dx , Dy in terms of the transverse quantities in Eqs. (1)–(13), we obtain a set of first-order non-homogeneous linear differential equations ðkÞ ;z

f

ðkÞ ðkÞ

¼C f

þ

ðkÞ bðkÞ f9mn ðzÞ;

f f

ðkÞ ;z

ðkÞ

b

¼ ¼ ¼

n n n

ðkÞ ðkÞ ðkÞ f1mn ðzÞf2mn ðzÞ    f8mn ðzÞ

o

ðkÞ    b8

o

0 6 6 0 6 6 ðkÞ 6c 6 31 6 ðkÞ 6c 6 41 ¼6 6 cðkÞ 6 51 6 ðkÞ 6c 6 61 6 6 0 4

0

c13

ðkÞ

c14

ðkÞ

c15

0

0

0

c23

ðkÞ

c24

ðkÞ

0

c26

ðkÞ

0

c32

ðkÞ

0

0

0

0

c37

c42

ðkÞ

0

0

0

0

c47

c52

ðkÞ

0

0

0

0

c57

c62

ðkÞ

0

0

0

0

c67

0

0

0

c75

ðkÞ

c76

ðkÞ

0

0

ðkÞ c84

ðkÞ c85

ðkÞ c86

0

0

0

ðkÞ

ðkÞ ðkÞ ðkÞ ðkÞ

0

3

7 0 7 7 7 ðkÞ c38 7 7 7 ðkÞ 7 c48 7 7: ðkÞ c58 7 7 7 ðkÞ 7 c68 7 7 0 7 5 0

The basis of solution space of the homogeneous form of Eq. (60) is

;

ðkÞ ðkÞ ðkÞ f1mn;z ðzÞf2mn;z ðzÞ    f8mn;z ðzÞ

ðkÞ ðkÞ b1 b2

C ðkÞ

ð60Þ

where ðkÞ

2

o

ðkÞ

ðkÞ

ðkÞ

f i ¼ Di eqi z ;

ð61Þ

; where i ¼ 1, 2, . . ., 8, and qimn are the characteristic roots ðkÞ of the matrix C ðkÞ ; Di are the corresponding characteristic vectors.

1208

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

Fig. 7. rz through thickness of the plate in loading: (a) case 1 (x ¼ a=2, y ¼ b=2), (b) case 2 (x ¼ a=2, y ¼ b=2) and (c) case 3 (x ¼ a=2, y ¼ b=2).

The roots come in pairs with the same absolute value and there are usually two pairs of conjugated complex roots and two pairs of real roots for piezoelectric material and three-pair real roots for common fiber-reinforced composite material. If the characteristic roots are complex number, in order to get the solution in real variables, the components of solution space correðkÞ sponding to each pair of complex roots, qimn ¼ aimn þ ðkÞ ibimn , can be derived as 
 h
i ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ Re Di cos bimn z
Im Di sin bimn z eaimn z ; ð62Þ h



 i ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ Im Di cos bimn z þ Re Di sin bimn z eaimn z ;

ð63Þ

where, Re and Im denote the real and imaginary parts of a complex quantity. In processing general solution, there are eight undetermined unknowns for each of eight components of basic solution space and the specified solution for Eq. (60) is f

ðkÞ

ðkÞ

¼
A




ðkÞ


1

where elements of the matrixes and vectors are given in appendix and I is 8  8 unit matrix. Specifically for the non-piezoelectric material layer, the equations containing the terms DðkÞ and ðkÞ in Eqs. z (1)–(54) may be removed and all electric permittivity coefficients eij and electric field as well as dielectric components are taken as zero; the corresponding matrix and coefficients are listed also in Appendix A. Finally there will be 6N þ 4 unknowns following with 6N þ 4 equations made of 6ðN
1Þ þ 2 interfaces conditions and 8 surface conditions. Therefore the solution can be obtained arithmetically. Notice that Eqs. (62)–(64) have a concise form and different from Pagano’s approach [1,5] in which a eighth-order algebraic equation must be solved and more sophisticated mathematical derivation must be made.

4. Results and discussion ðkÞ bðkÞ f9mn ðzÞ

cmn I þ C
ðkÞ 
1 ðkÞ þ BðkÞ cmn I
C bðkÞ f9mn ðzÞ;

ð64Þ

A specific laminate plate presented here consists of five physical layers, laid in fiber-orientation 0°/0°/90°/0°/ 0° with the top and bottom, PZT-5A piezoelectric ma-

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

1209

Fig. 8. In-plane shear stress, rx , ry through thickness of the plate in loading: (a) case 1 (x ¼ a=2, y ¼ b=2), (b) case 2 (x ¼ a=2, y ¼ b=2) and (c) case 3 (x ¼ a=2, y ¼ b=2).

terial layers and the three substrates, graphite-epoxy material layers. Both the length and breadth of the plate are 400 mm. The thickness of piezoelectric layer, h1 is 1.0 mm and the same thickness of three substrates, h2 alters from 4 to 40 mm in order to reflect the effect of the ratio of thickness/span and the ratio of h1 =h2 C11 ¼ C22 ¼ 99:201 ðGpaÞ;

C12 ¼ 54:016;

C13 ¼ C23 ¼ 50:778 ðGpaÞ;

C33 ¼ 86:856 ðGpaÞ;

C44 ¼ C55 ¼ 21:1 ðGpaÞ; e15 ¼ e24 ¼ 12:322 ðm2 =CÞ;

C66 ¼ 22:6 ðGpaÞ; e31 ¼ e32 ¼ 7:209 ðm2 =CÞ;

The material properties of the graphite–epoxy layers are in the process of computation, one term, namely n ¼ m ¼ 1, is truncated in the trigonometric series. The results reported in this paper are the mechanical responses in for three loading cases. Case (1) pure mechanical loading,
1z11 ¼ 1:0; q
q1xz11


2z11 ¼ 0:0; q

¼
q1yz11 ¼
q2yz11 ¼ 0:0; 1 2 /
111 ¼ /
211 ¼ 0:0; T
11 ¼ T
11 ¼ 0:0: ¼


q2xz11

e33 ¼ 15:1184 ðm2 =CÞ; e11 ¼ e22 ¼ 1:53ð10
8 farads=mÞ;

Case (2) mechanical and electric loading,

e33 ¼ 1:50ð10
8 farads=mÞ:


q1z11 ¼ 1:0;

The material properties of the piezoelectric layers are C11 ¼ 183:44 ðGpaÞ; C12 ¼ C13 ¼ 4:363 ðGpaÞ C22 ¼ C33 ¼ 11:662 ðGpaÞ; C23 ¼ 3:918 ðGpaÞ C44 ¼ 2:87 ðGpaÞ;

C55 ¼ C66 ¼ 7:17 ðGpaÞ


q2z11 ¼ 0:0;


q1xz11 ¼
q2xz11 ¼
q1yz11 ¼
q2yz11 ¼ 0:0; 1 2 1 2 /
11 ¼ 1:0; /
11 ¼ 0:0; T
11 ¼ T
11 ¼ 0:0: Case (3) mechanical, electric loading and temperature change,

1210

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212


q1z11 ¼ 1:0;
q1xz11


q2xz11


q2z11 ¼ 0:0;
q1yz11


q2yz11

¼ ¼ ¼ 1 2

/11 ¼ 1:0; /11 ¼ 0:0;

¼ 0:0;

1
pz11 ¼ 1:0;

2
pz11 ¼ 0:0:

The charts in Figs. 2–7 demonstrate the shear stresses, axial displacements, transverse normal stress and deflection of the plate subjected to three different loading cases. The scale of the ordinate of the piezoelectric layers are magnified to fit the charts and the solid line corresponds to H =L ¼ 3=10, the dash-dot line to H =L ¼ 3=25 and dot line to H =L ¼ 3=100. Notice that the graphics is drawn in spline line by using Matlab and the distribution charts of in-plane stresses in Fig. 8 reflect the discontinuous derivatives of the displacement solution. Tables 1–3 present the data of the shear stresses sxz and the axial displacement, u, of the plate with different H =L ratio in three different loading cases. The variations of transverse shear stresses are significant in the Limit State Analysis of the thin laminated

plate and Figs. 2 and 3 describe the distribution of transverse shear stresses sxz , syz , through the thickness of the plate. The shear stresses increase sharply in the substrate and peak in the vicinity of the mid-surface. The distribution of shear stress sxz induced by the applied potential is contrary to that by the mechanical loading. In contrast of shear stresses sxz and syz in Figs. 2 and 3b, it is observed that the shear response is sensitive to the material properties especially in the potential case in that there is a different fiber direction between the second, third, and fourth substrate layers of the plate. Figs. 2 and 3c demonstrate more important phenomena is that the that shear stress corresponding to the temperature change distribution oscillates through the thickness of the plate when the plate thickens in the combined loading. It is observed from Figs. 4 and 5 that the plate mainly undergoes the linear bending deformation in extreme thinness. The induced deformation by the applied potential and the temperature change increases and is nonlinear when the plate becomes thicker.

Table 1 Axial deformation and shear stress distribution across the thickness at the edge of plate (x ¼ 0, y ¼ L=2) with h1 ¼ 0:001, h2 ¼ 0:004 and H =L ¼ 3=100 Height Z )0.0070 )0.0065 )0.0060 )0.0040 )0.0020 0.0000 0.0020 0.0040 0.0060 0.0065 0.0070

u

sxz Case 1

Case 2

Case 3

Case 1

Case 2

Case 3

)0.72740E)13 )0.59275E)10 )0.11887E)09 7.7864 12.393 12.605 12.388 7.7804 )0.73755E)02 )0.36831E)02 0.73275E)14

0.41237E)10 0.14497E)07 0.29047E)07 )1921.5 )3067.9 )3112.6 )3052.8 )1954.0 )79.963 )40.442 0.18190E)11

0.45475E)11 )12.410 )24.941 )1471.4 )2327.8 )2376.4 )2359.2 )1505.1 )59.859 )30.332 0.45475E)12

0.38797E)09 0.38493E)09 0.38199E)09 0.25167E)09 0.12303E)09 0.10062E)12 )0.12283E)09 )0.25148E)09 )0.38181E)09 )0.38186E)09 )0.38191E)09

)0.95951E)07 )0.95250E)07 )0.94552E)07 )0.62782E)07 )0.31430E)07 )0.14899E)08 0.28455E)07 0.59804E)07 0.91563E)07 0.90924E)07 0.90284E)07

)0.73170E)07 )0.72647E)07 )0.72140E)07 )0.47789E)07 )0.23748E)07 )0.77145E)09 0.22204E)07 0.46264E)07 0.70658E)07 0.69995E)07 0.69332E)07

Table 2 Axial deformation and shear stress distribution across the thickness at the edge of plate (x ¼ 0, y ¼ L=2) with h1 ¼ 0:001, h2 ¼ 0:016 and H =L ¼ 3=25 u

Height Z

sxz Case 1

Case 2

Case 3

Case 1

Case 2

Case 3

)0.0250 )0.0245 )0.0240 )0.0160 )0.0080 0.0000 0.0080 0.0160 0.0240 0.0245 0.0250

)0.16614E)14 )0.12035E)11 )0.24130E)11 1.9434 2.8758 2.9168 2.8601 1.9366 )0.73755E)02 )0.36831E)02 )0.11102E)15

)0.12225E)12 0.35815E)09 0.71864E)09 )586.55 )877.49 )883.30 )860.22 )619.37 )79.963 )40.442 0.18190E)11

0.11369E)11 )9.6092 )19.312 )582.41 )821.86 )1047.7 )1416.4 )979.24 )71.179 )35.998 )0.20464E)11

0.25390E)10 0.25216E)10 0.25043E)10 0.13980E)10 0.44520E)11 0.20073E)12 )0.40945E)11 )0.13700E)10 )0.24882E)10 )0.24769E)10 )0.24657E)10

)0.77003E)08 )0.76575E)08 )0.76153E)08 )0.43720E)08 )0.16000E)08 )0.43572E)09 0.75648E)09 0.35303E)08 0.67298E)08 0.60670E)08 0.54046E)08

)0.94958E)08 )0.96029E)08 )0.97110E)08 )0.57709E)08 )0.21337E)08 0.72388E)10 0.19438E)08 0.65120E)08 0.12709E)07 0.11813E)07 0.10918E)07

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

1211

Table 3 Axial deformation and shear stress distribution across the thickness at the edge of plate (x ¼ 0, y ¼ L=2) with h1 ¼ 0:001, h2 ¼ 0:04 and H=L ¼ 3=10 Height Z )0.0610 )0.0605 )0.0600 )0.0400 )0.0200 0.0000 0.0200 0.0400 0.0600 0.0605 0.0610

u

sxz Case 1

Case 2

Case 3

Case 1

Case 2

Case 3

)0.68348E)16 )0.15465E)12 )0.31081E)12 0.76543 0.87116 0.86929 0.84367 0.80480 )0.73755E)02 )0.36831E)02 )0.33307E)15

0.35927E)13 0.57648E)10 0.11564E)09 )285.03 )331.50 )327.55 )313.15 )332.51 )79.963 )40.442 12733E)10

)0.22737E)12 )3.6698 )7.3755 )97.280 )151.28 )753.99 )1895.7 )2051.9 )78.974 )39.900 )0.45475E)11

0.47802E)11 0.47321E)11 0.46842E)11 0.13783E)11 )0.86330E)12 0.45880E)12 0.15499E)11 )0.10692E)11 )0.49844E)11 )0.49345E)11 )0.48847E)11

)0.17830E)08 )0.17669E)08 )0.17509E)08 )0.54679E)09 0.24858E)09 )0.31696E)09 )0.77834E)09 0.86714E)10 0.12753E)08 0.62114E)09 )0.32693E)10

)0.18183E)08 )0.19611E)08 )0.21040E)08 )0.18507E)08 )0.10251E)08 0.76829E)09 )0.32633E)08 0.23385E)08 0.17386E)07 0.16372E)07 0.15358E)07

Fig. 6 illustrates the variation of transverse displacement across the thickness of the plate. The figure shows that the transverse displacement is not an approximately constant when the plate becomes thick. As for the transverse normal stress, it become more complex when the thermal loading is involved. Fig. 7 gives that the transverse normal stress is linear through the thickness of the plate in pure mechanical loading and the combination of mechanical and potential loading. In the thermal loading, the transverse normal stress varies drastically and non-linearly when the plate thickens. Fig. 8 gives the through-thickness distribution of inplane normal stress rx , ry . Under pure mechanical loading, in Fig. 8a, it is observed that the normal stresses are zero at the mid-plane and have antisymmetric distribution in the graphite–epoxy layers but the highly stretching stresses are induced on the piezoelectric layers. Fig. 8b indicates that the external potential induces a converse stress distribution compared with those by the mechanical loading. In the combined loading case 3, the distribution of in-plane normal stresses is drastically affected by the external thermal loading and varies non-linearly.

From the present exact solution, the evaluation on the ratio of thickness/span of the plate is a useful guide for developing satisfactorily approximate solution of piezoelectric laminated plates with complex boundaries. Acknowledgements This work is supported by ARC grant, A89800361. The first author acknowledges useful discussions with Professor Paul Heyliger at Colorado State University. Appendix A Definition of matrix elements in Eq. (60)
 ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ R1 ¼ e31 am
am C13 e33 =C33 ;
 ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ R2 ¼ bn C23 e33 =C33 þ e33 ; ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

2

ðkÞ

A new procedure for piezoelectric structures is directly developed from the basic differential equations by using the trigonometric series and is used to obtain exact solutions of thermoelectroelastic laminates with simply supported and grounded boundaries. Numerical results show that both the electric and the temperature change significantly affect the displacement distribution of the plate due to the interaction between the thermal, piezoelectric, and mechanical fields; in thermal loading, stresses vary drastically and non-linearly across the thickness of the plate when the plate becomes thicker.

ðkÞ

ðkÞ

R4 ¼ e32 bn
e33 C23 bn =e33 ; ðkÞ

ðkÞ

R5 ¼ e33
e33 =C33 ; ðkÞ

5. Conclusions

ðkÞ ðkÞ

R3 ¼ e31 am
am C13 e33 =e33 ;

ðkÞ ðkÞ

ðkÞ

ðkÞ

R6 ¼ C33 e33 =e33
e33 ;
. ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ R7 ¼ k3 e33 =e33 þ r3 R6 ;
. ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ R5 ; R8 ¼ k3 e33 =C33 þ r3 ðkÞ

ðkÞ ðkÞ ðkÞ c14 ¼
e15 am =C55 ; ðkÞ ðkÞ 1=C55 c23 ¼
bn ; ðkÞ ðkÞ ðkÞ ðkÞ
bn e24 =C44 ; c26 ¼ 1=C44 ; ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ
R3 =R6 ; c32 ¼
R4 =R6 ; ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ e33 =ðe31 R6 Þ; c38 ¼
1=R6 ; ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ
R1 =R5 ; c42 ¼
R2 =R5 ; ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ e33 =ðC33 R5 Þ; c48 ¼
1=R5 ;

c13 ¼
am ; ðkÞ c15 ðkÞ c24 ðkÞ c31 ðkÞ c37 ðkÞ c41 ðkÞ c47

¼ ¼ ¼ ¼ ¼ ¼

1212

C. Zhang et al. / Computers and Structures 80 (2002) 1201–1212

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

c51 ¼ C11 a2m þ e31 am R1 =R5 þ C13 am R3 =R6 þ C66 b2n ; ðkÞ c52 ðkÞ c57

¼

ðkÞ c58 ðkÞ c61 ðkÞ c62 ðkÞ c67 ðkÞ c68 ðkÞ c75

¼

ðkÞ c84 ðkÞ c85 ðkÞ b1 ðkÞ b3 ðkÞ b5

¼

ðkÞ

¼ ¼

ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ am bn C12 þ e31 am R2 =R5 þ C13 an R4 =R6 þ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ
e31 e33 am =ðC33 R5 Þ
C13 am e33 =ðe33 R6 Þ; ðkÞ ðkÞ ðkÞ ðkÞ e31 am =R5 þ C13 am =R6 ; ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ am bn C12 þ bn C23 R3 =R6 þ e32 bm R1 =R5 þ ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ C66 am bn ;

ðkÞ

ðkÞ ðkÞ

ðkÞ

ðkÞ

ðkÞ ðkÞ

ðkÞ

ðkÞ

ðkÞ

¼
bn C23 e33 =ðe33 R6 Þ
e32 e33 bm =ðC33 R5 Þ; ¼ ¼ ¼ ¼

ðkÞ ðkÞ ðkÞ c86 ¼ e24 bn =C44 ; ðkÞ b8 ¼ 0; ðkÞ R7 ;

ðkÞ ðkÞ

b6 ¼

ðkÞ


a2m e11
b2n e22 ;

ðkÞ ðkÞ e15 am =C55 ; ðkÞ ðkÞ b2 ¼ b7 ¼

¼ R6 b4 ¼ ¼

ðkÞ

ðk Þ ðkÞ ðkÞ k1 am
e31 am R7 ðk Þ ðkÞ ðkÞ k2 bn
C23 bn R6

ðkÞ

ðkÞ


C12 am R6 ; ðkÞ

ðkÞ


e32 bn R7

In the case of non-piezoelectric substrate material the above coefficients are reduced as ðkÞ

c14 ¼ 1=C55 ;

ðkÞ

c25 ¼ 1=C44 ;

c13 ¼
am ; c23 ¼
bn ; ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

c31 ¼ C13 am =C33 ; ðkÞ

ðkÞ

ðkÞ

c32 ¼ C13 bn =C33 ;

ðkÞ

c36 ¼ 1=C33 ;


2 ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ c41 ¼ C11 a2m
C13 a2m =C33 þ C66 am bn ; ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

c42 ¼
C13 C23 am bn =C33 þ C12 am bn þ C66 am bn ; ðkÞ

c46 ¼
C13 aðkÞ =C33 ; ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

c51 ¼ C12 am bn
C13 C23 am bn =C33 þ C66 am bn ;
2 ðkÞ ðkÞ ðkÞ ðkÞ ðkÞ c52 ¼ C22 b2n
C23 b2n =C33 þ C66 a2m ; ðkÞ

ðkÞ

ðkÞ ðkÞ

ðkÞ

c56 ¼
C23 bn =C33 c64 ¼ am c65 ¼ bn ; ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

b1 ¼ b2 ¼ b6 ¼ 0; ðkÞ

b3 ¼ k3 =C33 ;

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

ðkÞ

b5 ¼ k2 bn
C23 bn k3 =C33

References

C66 am bn ;

¼ b2n C22 þ bn C23 R4 =R6 þ e32 bm R2 =R5 þ C66 a2m ; ðkÞ ðkÞ ðkÞ ðkÞ bn C23 =R6 þ e32 bm =R5 ; ðkÞ am c76 ¼ bn ; ðkÞ ðkÞ ðkÞ ðkÞ
ðe33 am Þ2 =C55
ðe24 bn Þ2 =C44

ðkÞ

b4 ¼ k1 am
C13 am k3 =C33 ;

[1] Pagano NJ. Exact solutions for rectangular bidirectional composites and sandwich plates. J Comp Mater 1970;4: 20–4. [2] Ray MC, Rao KM, Samanta B. Exact analysis of coupled electroelastic behaviour of a piezoelectric plate under cylindrical bending. Comput Struct 1992;45(4):667–77. [3] Ray MC, Rao KM, Samanta B. Exact solution for static analysis of an intelligent structure under cylindrical bending. Comput Struct 1993;47(6):1031–42. [4] Ray MC, Bhattacharya R, Samanta B. Exact solutions for static analysis of intelligent structures. AIAA J 1993;31(9): 1684–91. [5] Heyliger P, Brooks S. Exact solutions for laminated piezoelectric plates in cylindrical bending. J Appl Mech 1996;63(4):903–10. [6] Heyliger P, Brooks S. Free vibration of piezoelectric laminates in cylindrical bending. Int J Solids Struct 1995;32(20):2945–60. [7] Senthil S, Batra RC. Cylindrical bending of laminated plates with distributed and segmented piezoelectric actuators/sensors. AIAA J 2000;38(5):857–67. [8] Xu KM, Noor AK, Tang Y. Three-dimensional solutions for coupled thermoelectroelastic response of multilayered plates. Comput Meth Appl Mech Engng 1995;126:355–71. [9] Dube GP, Kapuria S, Dumir PC. Exact piezothermoelastic solution of simply-supported orthotropic flat panel in cylindrical bending. Int J Mech Sci 1996;38(11):1161–77. [10] Sung KH, Charles K. Finite element analysis of composite structures containing distributed piezoceramic sensors and actuators. AIAA J 1993;30(3):771–80. [11] Hwang WS, Park HC. Finite element modeling of piezoelectric sensors and actuators. AIAA J 1993;31(5):930–7. [12] Chattopadhyay A, Li J, Gu H. Coupled thermo-piezoelectric-mechanical model for smart composite laminates. AIAA J 1999;37(12):1633–8. [13] Mackerle J. Smart materials and structures: FEM and BEM simulations. A bibliography (1997–1999). Finite Elem Anal Des 2001;37:71–83.