Three-dimensional analytical solution for functionally graded magneto–electro-elastic circular plates subjected to uniform load

Available online at www.sciencedirect.com

Composite Structures 83 (2008) 381–390 www.elsevier.com/locate/compstruct

Three-dimensional analytical solution for functionally graded magneto–electro-elastic circular plates subjected to uniform load X.Y. Li a, H.J. Ding a, W.Q. Chen a

a,b,*

Department of Civil Engineering, Zhejiang University, Hangzhou 310027, PR China b State Key Lab of CAD&CG, Zhejiang University, Hangzhou 310027, PR China Available online 25 May 2007

Abstract The problem of a functionally graded, transversely isotropic, magneto–electro-elastic circular plate acted on by a uniform load is considered. The displacements and electric potential are represented by appropriate polynomials in the radial coordinate, of which the coefficients depends on the thickness coordinate, and are called the generalized displacement functions. The governing equations as well as the boundary conditions for these generalized displacement functions are derived from the original equations of equilibrium for axisymmetric problems and the boundary conditions on the upper and lower surfaces of the plate. Explicit expressions are then obtained through a step-by-step integration scheme, with five integral constants determinable from the boundary conditions at the cylindrical surface in the Saint Venant’s sense. The analytical solution is suited to arbitrary variations of material properties along the thickness direction, and can be readily degenerated into those for homogeneous plates. A particular circular plate, with some material constants being the exponential functions of the thickness coordinate, is finally considered for illustration.  2007 Elsevier Ltd. All rights reserved. Keywords: Transversely isotropic; Functionally graded; Magneto–electro-elastic; Circular plates; Direct displacement method

1. Introduction Magneto–electro-elastic materials (MEEMs), as a special sort of smart materials, have received more and more attention in recent years [1–3]. Since MEEM possesses piezo-electric, piezo-magnetic and magneto-electric effects, thereby making the composite sensitive to elastic, electric and magnetic fields [4,5], the magnetic, electric and mechanical energies in it can be converted from one form to another [6–11]. Furthermore magneto–electro-elastic materials exhibit stronger coupling effects in comparison with the single phase piezo-electric or piezo-magnetic materials [12]. When the material properties vary continuously in space on the macroscopic scale, the materials are referred to as functionally graded materials (FGMs) [13]. Functionally graded * Corresponding author. Address: Department of Civil Engineering, Zhejiang University, Hangzhou 310027, PR China. Tel.: +86 571 87952267; fax: +86 571 87952165. E-mail address: [email protected] (W.Q. Chen).

0263-8223/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2007.05.006

magneto–electro-elastic materials (FGMEEMs) combine the characters of both MEEMs and FGMs. Consequently, studies related to the performance of FGMEEMs under the action of external loads are very popular at present. Li [14] investigated the problems of multi-inclusion and inhomogeneity in MEEMs. For transversely isotropic magneto–electro–thermo-elastic materials, Chen and Lee [15] developed a new type state space method, which was then applied to analyzing the bending of inhomogeneous (or functionally graded) plates. Chen et al. [16] examined the free vibration of functionally graded magneto–electro– thermo-elastic plates and found two separated classes of vibrations. The exact static solutions were derived by Pan and Han [17] through pseudo-Stroh formula and the propagator matrix method for a multi-layered FGMEEM plate with material properties in each sub-layer varying in the same exponential law along the thickness direction. Recently, Bhangale and Ganesan [18,19], employed series solutions in the plane of the plate along with the finite element method involving the generalized displacements as

382

X.Y. Li et al. / Composite Structures 83 (2008) 381–390

basic freedoms in the thickness direction to study the static and free vibration problems of FGMEEM rectangular plates. However, to the authors’ knowledge, no research on FGMEEM circular plates on the basis of exact threedimensional equations has been reported yet. Hence, this paper intends to develop a three-dimensional (3D) analysis of the elastic, electric and magnetic fields in transversely isotropic FGMEEM circular plates subjected to a uniform load. The direct displacement method is employed, which expresses the two components of the mechanical displacement, the electric potential and the magnetic potential in terms of appropriate polynomials in the radial coordinate. The coefficients of the polynomials are functions of the thickness coordinate variable to be determined. In this way, the original problem is transformed into a system of ordinary differential equations of the generalized displacement functions, i.e. the unknown coefficients in the polynomials. A step-by-step integration is executed carefully to find the explicit expressions of the generalized displacement functions, with the resultant integral constants determined from the boundary conditions at the cylindrical surface of the plate. The present solution satisfies exactly the boundary conditions on the upper and bottom surfaces of the plate, while approximately the boundary conditions on the cylindrical surface in the Saint Venant’s sense. The significant merit of the solution lies in the fact that the material parameters can vary in an arbitrary manner provided that the positive definiteness of the generalized strain energy and some integral conditions during the derivation are fulfilled. In this sense, the solution can be readily degenerated into that for homogeneous plates. A particular FGMEEM plate is considered and numerical results are presented to illustrate the effect of material inhomogeneity on the elastic, electric and magnetic fields in the plate. 2. Basic equations In a cylindrical coordinate system (r, h, z), the basic equations for axisymmetric problems of MEEMs [20] are listed as follows: rr;r þ srz;z þ r1 ðrr  rh Þ ¼ 0; srz;r þ r1 srz þ rz;z ¼ 0:

ð1Þ

Dr;r þ r1 Dr þ Dz;z ¼ 0:

ð2Þ

Br;r þ r1 Br þ Bz;z ¼ 0:

ð3Þ

rr ¼ c11 u;r þ c12 r1 u þ c13 w;z þ e31 /;z þ d 31 w;z ;

where the comma denotes differentiation with respect to the indicated variable; rr, rh, rz and srz are the stress components; Dr and Dz are the electric displacement components; Br and Bz are the magnetic induction components; u and w are the displacements in radial and thickness directions, respectively; / and w are the electric and magnetic potentials, respectively; cij, eij, eij, gij, lij and dij are, respectively, the elastic, piezo-electric, dielectric magneto-electric, magnetic and piezo-magnetic coefficients. In the following analysis, we assume that the material coefficients in Eqs. (4)–(6) are arbitrary functions of the thickness coordinate; for homogeneous materials, these coefficients are simply constant. It is noted that the magnetic, electric and elastic fields can be decoupled by setting the appropriate material coefficients zero in Eqs. (4)–(6). For example, the magnetic (electric) field becomes independent of the electro-elastic (magneto-elastic) coupling field if dij = 0 (eij = 0) and gij = 0. In particular, all the three fields are decoupled from each other when eij, dij and gij vanish identically. 3. Uniformly loaded circular plate Consider a circular plate with radius a and thickness h, subject to a uniform load q, as shown in Fig. 1. We will adopt the direct displacement method to solve this axisymmetric problem. In view of the fact that the radial displacement u vanishes at r = 0, we take u ¼ ru1 ðzÞ þ r3 u3 ðzÞ; w ¼ w0 ðzÞ þ r2 w2 ðzÞ þ r4 w4 ðzÞ;

ð7Þ

2

ð8Þ

2

ð9Þ

/ ¼ /0 ðzÞ þ r /2 ðzÞ; w ¼ w0 ðzÞ þ r w2 ðzÞ;

where ui(z) and wi(z) are referred to as the displacement functions, /i(z) as the electric potential functions and wi(z) as the magnetic potential functions; or the three are also holistically called as the generalized functions. The form assumed in Eqs. (7)–(9) extends that for homogeneous piezo-electric and MEEM materials [20,21]. Substituting Eqs. (7)–(9) into Eqs. (4)–(6) yields rr ¼ ðc11 þ c12 Þu1 þ c13 w0;z þ e31 /0;z þ d 31 w0;z þ r2 ½ð3c11 þ c12 Þu3 þ c13 w2;z þ e31 /2;z þ d 31 w2;z  þ r4 c13 w4;z ; rh ¼ ðc11 þ c12 Þu1 þ c13 w0;z þ e31 /0;z þ d 31 w0;z þ r2 ½ðc11 þ 3c12 Þu3 þ c13 w2;z þ e31 /2;z þ d 31 w2;z  þ r4 c13 w4;z ;

rh ¼ c12 u;r þ c11 r1 u þ c13 w;z þ e31 /;z þ d 31 w;z ; rz ¼ c13 ðu;r þ r1 uÞ þ c33 w;z þ e33 /;z þ d 33 w;z ; szr ¼ c44 ðu;z þ w;r Þ þ e15 /;r þ d 15 w;r :

ð4Þ

Dz ¼ e31 ðu;r þ r1 uÞ þ e33 w;z  e33 /;z  g33 w;z :

h /2

ð5Þ

Br ¼ d 15 ðu;z þ w;r Þ  g11 /;r  l11 w;r ; Bz ¼ d 31 ðu;r þ r1 uÞ þ d 33 w;z  g33 /;z  l33 w;z ;

h /2

o

Dr ¼ e15 ðu;z þ w;r Þ  e11 /;r  g11 w;r ;

a z

ð6Þ

Fig. 1. Circular plate subject to a uniform load.

r

X.Y. Li et al. / Composite Structures 83 (2008) 381–390

rz ¼ 2c13 u1 þ c33 w0;z þ e33 /0;z þ d 33 w0;z 2

4

þ r ½4c13 u3 þ c33 w2;z þ e33 /2;z þ d 33 w2;z  þ r c33 w4;z ;

383

½c44 ðu3;z þ 4w4 Þz¼h=2 ¼ 0;

ð28Þ

½c44 ðu1;z þ 2w2 Þ þ 2e15 /2 þ 2d 15 w2 z¼h=2 ¼ 0;

ð29Þ

½4e31 u3 þ e33 w2;z  e33 /2;z  g33 w2;z z¼h=2 ¼ 0;

ð30Þ

½2e31 u1 þ e33 w0;z  e33 /0;z  g33 w0;z z¼h=2 ¼ 0;

ð31Þ

½4d 31 u3 þ d 33 w2;z  g33 /2;z  l33 w2;z z¼h=2 ¼ 0;

ð32Þ

½2d 31 u1 þ d 33 w0;z  g33 /0;z  l33 w0;z z¼h=2 ¼ 0:

ð33Þ

3

srz ¼ r½c44 ðu1;z þ 2w2 Þ þ 2e15 /2 þ 2d 15 w2  þ r c44 ðu3;z þ 4w4 Þ: ð10Þ Dr ¼ r½e15 ðu1;z þ 2w2 Þ  2e11 /2  2g11 w2  þ r3 e15 ðu3;z þ 4w4 Þ; Dz ¼ 2e31 u1 þ e33 w0;z  e33 /0;z  g33 w0;z 2

4

þ r ð4e31 u3 þ e33 w2;z  e33 /2;z  g33 w2;z Þ þ r e33 w4;z : ð11Þ 3

Br ¼ r½d 15 ðu1;z þ 2w2 Þ  2g11 /2  2l11 w2  þ r d 15 ðu3;z þ 4w4 Þ; Bz ¼ 2d 31 u1 þ d 33 w0;z  g33 /0;z  l33 w0;z þ r2 ð4d 31 u3 þ d 33 w2;z  g33 /2;z  l33 w2;z Þ þ r4 d 33 w4;z : ð12Þ Introducing Eqs. (10)–(12) into Eqs. (1)–(3) and letting the coefficients of the each power of r vanish, we arrive at the following 11 differential equations governing the generalized displacement functions ðc33 w4;z Þ;z ¼ 0; ðe33 w4;z Þ;z ¼ 0; ðd 33 w4;z Þ;z ¼ 0;

ð13Þ

½c44 ðu3;z þ 4w4 Þ;z þ 4c13 w4;z ¼ 0;

ð14Þ

½c44 ðu1;z þ 2w2 Þ þ 2e15 /2 þ 2d 15 w2 ;z þ 8c11 u3 þ 2ðc13 w2;z þ e31 /2;z þ d 31 w2;z Þ ¼ 0; ð15Þ ð4c13 u3 þ c33 w2;z þ e33 /2;z þ d 33 w2;z Þ;z þ 4c44 ðu3;z þ 4w4 Þ ¼ 0; ð16Þ ð2c13 u1 þ c33 w0;z þ e33 /0;z þ d 33 w0;z Þ;z þ 2½c44 ðu1;z þ 2w2 Þ þ 2e15 /2 þ 2d 15 w2  ¼ 0: ð17Þ ð4e31 u3 þ e33 w2;z  e33 /2;z  g33 w2;z Þ;z þ 4e15 ðu3;z þ 4w4 Þ ¼ 0; ð18Þ ð2e31 u1 þ e33 w0;z  e33 /0;z  g33 w0;z Þ;z þ 2½e15 ðu1;z þ 2w2 Þ  2e11 /2  2g11 w2  ¼ 0; ð19Þ ð4d 31 u3 þ d 33 w2;z  g33 /2;z  l33 w2;z Þ;z þ 4d 15 ðu3;z þ 4w4 Þ ¼ 0; ð20Þ ð2d 31 u1 þ d 33 w0;z  g33 /0;z  l33 w0;z Þ;z þ 2½d 15 ðu1;z þ 2w2 Þ  2g11 /2  2l11 w2  ¼ 0;

ð21Þ

The boundary conditions of the circular plate at z = ±h/2 are z ¼ h=2 : rz ¼ q; z ¼ h=2 : rz ¼ 0;

srz ¼ 0;

srz ¼ 0;

Dz ¼ 0;

Dz ¼ 0;

Bz ¼ 0:

Bz ¼ 0:

4. Determination of the generalized displacement functions At the first sight, it seems difficult to determine the nine generalized functions in explicit forms, since the ordinary differential equations in Eqs. (13)–(21) are coupled and have variable coefficients. However, one can find that the boundary conditions in Eqs. (24)–(33) correspond exactly to certain terms in Eqs. (13)–(21). This amazing feature enables us to solve the system of differential equations step-by-step as described below. Integrating the three differential equations in Eq. (13) once from the lower limit z = h/2 and using the corresponding boundary conditions in Eq. (27), one obtains w4 ¼ a1 ;

ð34Þ

where a1 is an integral constant. The expression for u3 is obtained by integrating Eq. (14) and utilizing Eq. (34) and the boundary condition in Eq. (28) u3 ¼ 4a1 z þ a2 ;

ð35Þ

where a2 is another integral constant. Integrating Eqs. (16), (18) and (20) and employing Eqs. (34) and (35) and the boundary conditions in Eqs. (26), (30) and (32) leads to simultaneous equations of w2,z, /2,z and w2,z. Further integration leads to w2 ¼ 16a1 f21 ðzÞ  4a2 f20 ðzÞ þ a3 ;

ð36Þ

/2 ¼ 16a1 g21 ðzÞ  4a2 g20 ðzÞ þ a4 ;

ð37Þ

w2 ¼ 16a1 h21 ðzÞ  4a2 h20 ðzÞ þ a5 ;

ð38Þ

where a3, a4 and a5 are three new integral constants and Z z ½f2j ðzÞ; g2j ðzÞ; h2j ðzÞ ¼ nj ½f2 ðnÞ; g2 ðnÞ; h2 ðnÞ dn; h=2

ð22Þ

ðj ¼ 0; 1Þ;

ð23Þ T

Substitution of the expressions for rz, srz, Dz and Bz in Eqs. (10)–(12) into the Eqs. (22) and (23) results in the boundary conditions compatible with the system of differential equations in Eqs. (13)–(21) as follows:

ð39Þ T

½f2 ðzÞ; g2 ðzÞ; h2 ðzÞ ¼ C½c13 ; e31 ; d 31  ; 2 3 l33 e33  g233 e33 l33  d 33 g33 d 33 e33  e33 g33 16 7 C ¼ 4 l33 e33  d 33 g33 ðc33 l33 þ d 233 Þ c33 g33 þ e33 d 33 5; J d 33 e33  e33 g33 c33 g33 þ e33 d 33 ðc33 e33 þ e233 Þ

½2c13 u1 þ c33 w0;z þ e33 /0;z þ d 33 w0;z z¼h=2 ¼ q;

ð24Þ

½2c13 u1 þ c33 w0;z þ e33 /0;z þ d 33 w0;z z¼h=2 ¼ 0;

ð25Þ

½4c13 u3 þ c33 w2;z þ e33 /2;z þ d 33 w2;z z¼h=2 ¼ 0;

ð26Þ

J ¼ ðc33 l33 þ

ð27Þ

Introducing Eqs. (35)–(38) into Eq. (15), integrating once and taking advantage of the boundary condition in Eq. (29), one derives

½c33 w4;z z¼h=2 ¼ 0; ½e33 w4;z z¼h=2 ¼ 0; ½d 33 w4;z z¼h=2 ¼ 0;

ð40Þ d 233 Þe33



c33 g233

þ

l33 e233

 2e33 g33 d 33 :

384

X.Y. Li et al. / Composite Structures 83 (2008) 381–390

c44 ðu1;z þ 2w2 Þ þ 2e15 /2 þ 2d 15 w2 ¼ 32a1 H 1 ðzÞ  8a2 H 0 ðzÞ;

ð41Þ

4a1 H 1 ðh=2Þ  a2 H 0 ðh=2Þ ¼ 0;

ð42Þ

where H i ðzÞ ¼

Z

in which

ni ½c11  c13 f2 ðnÞ  e31 g2 ðnÞ  d 31 h2 ðnÞ dn;

h=2

u1 ¼ 32a1 F 1 ðzÞ  8a2 F 2 ðzÞ  2a3 z  2K 0 a4  2K 1 a5 þ a6 ; ð44Þ where a6 is another integral constant and

Z

h=2 z

h=2 z

ðe215 c1 44 þ e11 Þ dn; T 12 ðzÞ ðd 15 e15 c1 44 þ g 11 Þ dn;

ð53Þ

½K 0 ðnÞH i ðnÞ  T 11;n ðnÞg2i ðnÞ

 T 12;n ðnÞh2i ðnÞ dn;

ði ¼ 0; 1Þ:

Substituting Eqs. (36)–(38) and (44) into Eq. (21), integrating once and using the boundary condition in Eq. (33) results in d 33 w0;z  g33 /0;z  l33 w0;z ¼ 64a1 ½E1 ðzÞ þ d 31 F 1 ðzÞ þ 16a2 ½E0 ðzÞ þ d 31 F 0 ðzÞ

K 1 ¼ d 15 c1 44 ;

F i ðzÞ ¼ c1 44 H i ðzÞ  f2i ðzÞ  K 0 g2i ðzÞ  K 1 h2i ðzÞ; Z z ½F i ðnÞ; K i ðnÞ dn; ði ¼ 0; 1Þ: ½F i ðzÞ; K i ðzÞ ¼

z

h=2

ð43Þ

Substituting Eqs. (36)–(38) into Eq. (41) and integrating once, gives

K 0 ¼ e15 c1 44 ;

Z

¼ Z Gi ðzÞ ¼

z

ði ¼ 0; 1Þ:

T 11 ðzÞ ¼

ð45Þ

þ 4a3 d 31 z þ 4a4 ½T 12 ðzÞ þ d 31 K 0 ðzÞ þ 4a5 ½T 22 ðzÞ þ d 31 K 1 ðzÞ  2d 31 a6 :

ð55Þ

h=2

Introducing Eq. (41) into Eq. (17) and integrating once leads to 2c13 u1 þ c33 w0;z þ e33 /0;z þ d 33 w0;z ¼ 64a1 H 1 ðzÞ þ 16a2 H 0 ðzÞ þ a7

ð46Þ

ð47Þ

Substitution of Eq. (46) into the boundary conditions in Eqs. (24) and (25) yields a7 ¼ q;

ð48Þ

4a1 H 1 ðh=2Þ  a2 H 0 ðh=2Þ ¼ q=16:

ð49Þ

ð56Þ

ði ¼ 0; 1Þ:

Solving Eqs. (50), (51) and (54) simultaneously yields the expressions for w0,z, /0,z and w0,z, integration of which yields 5 X ai P 0i ðzÞ w0 ¼ 64a1 P 01 ðzÞ þ 16a2 P 02 ðzÞ þ 4 i¼3

 2a6 P 06 ðzÞ  qP 07 ðzÞ þ a8 ;

ð57Þ

/0 ¼ 64a1 Q01 ðzÞ þ 16a2 Q02 ðzÞ þ 4

Introducing Eq. (44) into Eq. (46) gives arise to

5 X

ai Q0i ðzÞ

i¼3

 2a6 Q06 ðzÞ  qQ07 ðzÞ þ a9 ;

c33 w0;z þ e33 /0;z þ d 33 w0;z ¼  64a1 ½H 1 ðzÞ þ c13 F 1 ðzÞ 16a2 ½H 0 ðzÞ

þ c13 F 0 ðzÞ 4a4 c13 K 0 ðzÞ

ð58Þ

w0 ¼ 64a1 R01 ðzÞ þ 16a2 R02 ðzÞ þ 4

5 X

ai R0i ðzÞ

i¼3

þ 4a5 c13 K 1 ðzÞ  2c13 a6  q: ð50Þ Substituting Eqs. (36)–(38) and (44) into Eq. (19), integrating once and making use of the boundary condition in Eq. (31), we get

 2a6 R06 ðzÞ  qR07 ðzÞ þ a0 ;

ð59Þ

where a0, a8 and a9 are integral constants, and Z z ½P 0i ðzÞ; Q0i ðzÞ; R0i ðzÞ ¼ ½P i ðnÞ; Qi ðnÞ; Ri ðnÞ dn; h=2

ði ¼ 1; 2; . . . ; 7Þ: T

Vi ¼ CSi ; Vi ¼ ½P i ðzÞ; Qi ðzÞ; Ri ðzÞ ; ði ¼ 1; 2; . . . ; 7Þ;

e33 w0;z  e33 /0;z  g33 w0;z

S1 ¼ ½H 1 ðzÞ þ c13 F 1 ðzÞ; G1 ðzÞ þ e31 F 1 ðzÞ; E1 ðzÞ þ d 31 F 1 ðzÞT ;

¼ 64a1 ½G1 ðzÞ þ e31 F 1 ðzÞ þ 16a2 ½G0 ðzÞ þ e31 F 0 ðzÞ

T

S2 ¼ ½H 0 ðzÞ þ c13 F 0 ðzÞ; G0 ðzÞ þ e31 F 0 ðzÞ; E1 ðzÞ þ d 31 F 0 ðzÞ ;

þ 4a3 e31 z þ 4a4 ½T 11 ðzÞ þ e31 K 0 ðzÞ þ 4a5 ½T 12 ðzÞ þ

Z z T 22 ðzÞ ¼ ðd 215 c1 44 þ l11 Þ dn; h=2 Z z Ei ðzÞ ¼ ½K 1 ðnÞH i ðnÞ  T 12;n ðnÞg2i ðnÞ  T 22;n ðnÞh2i ðnÞ dn;

h=2

e31 K 1 ðzÞ

where

h=2

where a7 is a integral constants, and Z z  H i ðzÞ ¼ H i ðnÞ dn; ði ¼ 0; 1Þ:

þ þ 4a3 c13 z þ

ð54Þ

a4 T 12 ðh=2Þ þ a5 T 22 ðh=2Þ ¼ 16a1 E1 ðh=2Þ  4a2 E0 ðh=2Þ;

 2e31 a6 :

ð51Þ

a4 T 11 ðh=2Þ þ a5 T 12 ðh=2Þ ¼ 16a1 G1 ðh=2Þ  4a2 G0 ðh=2Þ;

S4 ¼ ½c13 K 0 ðzÞ; e31 K 0 ðzÞ þ T 11 ðzÞ; d 31 K 0 ðzÞ þ T 12 ðzÞT ; T

S5 ¼ ½c13 K 1 ðzÞ; e31 K 1 ðzÞ þ T 12 ðzÞ; d 31 K 1 ðzÞ þ T 22 ðzÞ ; T

T

T

S3 ¼ ½c13 z; e31 z; d 31 z ; S6 ¼ ½c13 ; e31 ; d 31  ; S7 ¼ ½1; 0; 0 : ð52Þ

ð60Þ

X.Y. Li et al. / Composite Structures 83 (2008) 381–390

The above procedure yields explicit expressions for the generalized displacement functions, which involve a total of 10 integral constants ai ði ¼ 0; 1; 2;    ; 9Þ. Among them, a7 has been given by Eq. (48), a1 and a2 are given by Eqs. (42) and (49), and a4 and a5 by Eqs. (52) and (55). There are still five constants, i.e. a3, a6, a8, a9 and a0, to be determined, generally from the boundary conditions on the cylindrical surface (r = a). However, it can be shown that a8 is a rigid-body translation in the axial direction while a9 and a0 are pertinent to the reference electric and magnetic potentials, respectively. Therefore, these three constants do not alter the distribution of the stress, the electric displacement as well as the magnetic induction.

N 2j ¼

Z

h=2

h=2

N 3j

385

zj ½ðc11 þ c12 ÞF 0 ðzÞ  2c13 P 2 ðzÞ

 2e31 Q2 ðzÞ  2d 31 R2 ðzÞ dz; Z h=2 ¼ zj ½ðc11 þ c12 Þz  2c13 P 3 ðzÞ h=2

N 4j

 2e31 Q3 ðzÞ  2d 31 R3 ðzÞ dz; Z h=2 ¼ zj ½ðc11 þ c12 ÞK 0 ðzÞ  2c13 P 4 ðzÞ

N 5j

 2e31 Q4 ðzÞ  2d 31 R4 ðzÞ dz; Z h=2 ¼ zj ½ðc11 þ c12 ÞK 1 ðzÞ  2c13 P 5 ðzÞ

N 6j

 2e31 Q5 ðzÞ  2d 31 R5 ðzÞ dz; Z h=2 ¼ zj ½ðc11 þ c12 Þ  2c13 P 6 ðzÞ

h=2

ð64Þ

h=2

5. Determination of the integral constants Substituting back the generalized displacement functions into Eqs. (7)–(12) yields the expressions for displacements, electric and magnetic potentials, stresses, electric displacements and magnetic inductions for the problem, as listed in Appendix. It is seen that: (1) The expressions for the axial stress (rz), shear stress (srz), electric displacements (Dr, Dz) and magnetic inductions (Br, Bz) contain no undetermined constants, and hence they have been completely determined. This means that the circumferential boundary condition has no influence on these field variables; (2) rz, Dz and Bz are functions of the variable z only and independent of the radial coordinate r; (3) srz, Dr and Br vary linearly with r; (4) The radial and hoop stresses are equal to each other along the axisymmetric axis; (5) The average electric displacement Dr and magnetic induction Br are zero, i.e., Z Z 1 h=2 1 h=2   Dr  Dr dz ¼ 0; Br  Br dz ¼ 0; ð61Þ h h=2 h h=2 where Eqs. (52) and (55) have been used. The radial resultant force N(r) and bending moment M(r) are calculated as follows: Z h=2 rr dz ¼ 4a1 ½r2 L1  8N 10  þ a2 ½r2 L0  8N 20  N ðrÞ  h=2

ð62Þ  2a3 N 30  2a4 N 40  2a5 N 50 þ a6 N 60  qN 70 ; Z h=2 MðrÞ  zrr dz ¼ 4a1 ½r2 L2  8N 11  þ a2 ½r2 L1  8N 21  h=2

 2a3 N 31  2a4 N 41  2a5 N 51 þ a6 N 61  qN 71 ; where Z Li ¼

N 7j

 2e31 Q6 ðzÞ  2d 31 R6 ðzÞ dz; Z h=2 ¼ zj ½c13 P 7 ðzÞ þ e31 Q7 ðzÞ h=2

þ d 31 R7 ðzÞ dz;

h=2

zi ½ð3c11 þ c12 Þ  4c13 f2 ðzÞ  4e31 g2 ðzÞ

 4d 31 h2 ðzÞ dz; ði ¼ 0; 1; 2Þ; Z h=2 ¼ zj ½ðc11 þ c12 ÞF 1 ðzÞ  2c13 P 1 ðzÞ h=2

 2e31 Q1 ðzÞ  2d 31 R1 ðzÞ dz;

ðj ¼ 0; 1Þ:

For a simply-supported circular plate, the boundary conditions are N ðaÞ ¼ 0;

MðaÞ ¼ 0;

wða; 0Þ ¼ 0;

Dr ðaÞ ¼ 0;

Br ðaÞ ¼ 0:

ð65Þ

where the last two conditions have been satisfied as shown by Eq. (61). Introducing Eqs. (62) and (63) into the first two conditions in Eq. (65) yields 2a3 N 30  a6 N 60 ¼ 4a1 ½a2 L1  8N 10  þ a2 ½a2 L0  8N 20   2a4 N 40  2a5 N 50  qN 70 ; 2

ð66Þ

2

2a3 N 31  a6 N 61 ¼ 4a1 ½a L2  8N 11  þ a2 ½a L1  8N 21   2a4 N 41  2a5 N 51  qN 71 ; ð67Þ from which a3 and a6 can be determined. The constant a8 can be fixed by introducing Eq. (A10) into the third boundary condition in Eq. (65). Once the reference electric and magnetic potentials are prescribed, the electric and magnetic potentials will be completely determined. When the circular plate is clamped, the boundary conditions at r = a are uða; 0Þ ¼ 0;

wða; 0Þ ¼ 0;

 1; /ða; 0Þ ¼ /

h=2

N 1j

ð63Þ

h=2

w;r ða; 0Þ ¼ 0;  1: wða; 0Þ ¼ w

ð68Þ

In this case, a3 and a6 are determined by introducing Eqs. (A10) and (A9) into the third and first boundary conditions in Eq. (68), respectively, while the expressions for a8, a9 and a0 can be derived by substituting Eqs. (A10)–(A12) into the second, fourth and fifth boundary conditions. Substituting the integral constants such fixed back into Eqs. (A1)–(A12), we arrive at the 3D magneto–electroelastic analytical solutions for uniformly loaded FGMEEM

386

X.Y. Li et al. / Composite Structures 83 (2008) 381–390

circular plates with simply supported or clamped edges. These solutions bear the following merits:

Table 2 Comparison of deflection W(0,0) and moment Mð0Þ b = 0.01

(1) The material coefficients can be the arbitrary functions of z, and hence they can be applied to various FGM models or the homogeneous case. (2) They can be used to analyze uniformly loaded plates of piezo-magnetic, piezo-electric or elastic materials, if appropriate material coefficients are set zero. (3) The elastic, electric and magnetic fields are expressed explicitly, showing clearly the effect of material inhomogeneity and facilitating greatly the optimization design of FGM structures.

6. Numerical analysis In the numerical calculation, we take q = 1 · 106 (N/m2) and a = 0.1 (m). The following dimensionless quantities are also introduced for the convenience of display: r 1 z 1 rr 6 1;  6 g ¼ 6 ; rn ¼ ; a 2 h 2 q rz srz rg ¼ ; sng ¼ ; q q Dr Dz Dn ¼ pffiffiffiffiffiffiffiffiffiffiffi ; Dg ¼ pffiffiffiffiffiffiffiffiffiffiffi ; 0 0 c33 e33 c033 e033 Br Bz Bn ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ; Bg ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ; 0 0 c33 l33 c033 l033 sffiffiffiffiffiffi u w / e033 ; U ¼ ; W ¼ ; U¼ a h h c033 sffiffiffiffiffiffiffi w l033 h M ; b¼ ; M ¼ 2; W¼ 0 h c33 a qh

Ref. [22] Present

b = 0.1

W(0,0)

Mð0Þ

W(0,0)

Mð0Þ

6.2469 · 102 6.2421 · 102

2220 · 104 2216 · 104

6.2916 · 102 6.2916 · 102

2220 · 102 2216 · 102

Example 2. Now consider a particular circular FGMEEEM plate (h = 0.01) with material coefficients taken in the following form: c11 ¼ c011 ekðgþ1=2Þ ;

c12 ¼ c012 ekðgþ1=2Þ ;

c33 ¼ c033 ekðgþ1=2Þ ;

c44 ¼ c044 ;

e31 ¼ e031 ekðgþ1=2Þ ;

e33 ¼ e033 ekðgþ1=2Þ ;

d 31 ¼ d 031 ekðgþ1=2Þ ;

d 33 ¼ d 033 ekðgþ1=2Þ ;

d 15 ¼ d 015 ;

e11 ¼ e011 ;

c13 ¼ c013 ekðgþ1=2Þ ; e15 ¼ e015 ;

e33 ¼ e033 ekðgþ1=2Þ

g11 ¼ g011 ;

l33 ¼ l033 : ð70Þ

06n¼

c0ij ,

ð69Þ

We also assume that /(a,0) = 0 and w(a,0) = 0. Example 1. To validate the derivation as well as programming, we first consider a simply-supported homogeneous PZT-4 circular plate. The material constants of PZT-4 are listed in Table 1. Table 2 gives the dimensionless central deflection W(0,0) and central moment M(0) for b = 0.1 and 0.01. Our results, which are obtained based on the present analytical solutions by setting dij = gij = 0 and letting all other material coefficients constant, are compared to those given in Ref. [22]. We can see that a good agreement can be obtained, thus validating the current solutions.

e0ij ,

d 0ij ,

e0ij ,

The material property at z = h/2 is identical to CoFe2O4 [14,20], whose material constants are given in Table 3. As mentioned earlier, rg, sng, Dg, Dn, Bg and Bn are independent of the circumferential conditions, i.e. they are identical for simply-supported conditions and clamped conditions. Figs. 2, 4 and 6 depict the distributions of axial stress rg, electric displacement Dg and magnetic induction Bg, respectively, while those for sng, Dn and Bn at n = 1/2 are given in Figs. 3, 5 and 7, respectively. Fig. 2 shows that the absolute value of the normal stress increases with the material gradient index k. From Fig. 3,

Table 3 Material properties of CoFe2O4 2

Elastic (10 N/m ) PZT-4

Elastic (1010 N/m2)

c011 ¼ 12:6; c012 ¼ 7:78; c013 ¼ 7:43; c033 ¼ 11:5; c044 ¼ 2:56 e015 ¼ 12:7; e031 ¼ 5:2; e033 ¼ 15:1 e011 ¼ 646:3; e033 ¼ 562:2

Piezo-electric (C/m2) Dielectric (1011 F/m)

CoFe2O4 9

Property

l0ij

where and are the material constants at the surface z = h/2, g = z/h(1/2 6 g 6 1/2) is the dimensionless axial coordinate variable, and k is a dimensionless gradient index, which determines the degree of variation of material properties along the thickness. It is obvious that k = 0 corresponds to a homogeneous material. It is shown that the expressions for the elastic, electric and magnetic fields can be derived explicitly in terms of elementary functions; however, they are omitted here for the sake of simplicity.

Property Table 1 Material properties of PZT-4

g0ij

Dielectric (1011 F/m) Magnetic (106 Ns2/C2) Piezo-electric (C/m2) Magneto-electric (1012 Ns/VC) Piezo-magnetic (C/m2)

c011 ¼ 286:0; c012 ¼ 173:0; c013 ¼ 170:5; c033 ¼ 269:5; c044 ¼ 45:3 e011 ¼ 8:0; e033 ¼ 9:3 l011 ¼ 590; l033 ¼ 157 e015 ¼ 11:60; e031 ¼ 4:40; e033 ¼ 18:60 g011 ¼ 5:0; g033 ¼ 3:0 d 015 ¼ 550;

d 031 ¼ 580:3;

d 033 ¼ 699:7

X.Y. Li et al. / Composite Structures 83 (2008) 381–390

387

1.5

0 -0.1

1 -0.2 0.5

-0.3

Dξ (×10-4)

ση

-0.4 -0.5 -0.6

k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

-0.7 -0.8

-0.5 k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

-1 -1.5

-0.9 -1 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

-2 -0.5

0.5

η

-0.4

-0.1

0

0.1

0.2

0.3

0.4

0.5

1

-0.5

0

k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

-1

-1

Bη (×10-7)

-1.5 -2 -2.5

-2 -3 k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

-3 -4

-3.5 -5

-4 -4.5 -0.5 -0.4

-0.3

-0.2

-0.1

0

0.1

η

0.2

0.3

0.4

-6 -0.5

0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

η

Fig. 3. Dimensionless shear stress (n = 1/2).

Fig. 6. Dimensionless axial magnetic induction.

4

4 k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

2 0

2 0

Bξ (×10-6)

-2

Dη (×10-6)

-0.2

Fig. 5. Dimensionless radial electric displacement (n = 1/2).

0

-4 -6

-2 -4

k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

-6

-8

-8

-10 -12 -0.5

-0.3

η

Fig. 2. Dimensionless axial stress.

τξη

0

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

η Fig. 4. Dimensionless axial electric displacement.

0.4

0.5

-10 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

η Fig. 7. Dimensionless radial magnetic induction.

388

X.Y. Li et al. / Composite Structures 83 (2008) 381–390 400

4 k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

300 200 2

U (10-4)

σξ

100 0

-100 -200 -300 -400 -0.5

0

k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

-2 -0.5 -0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

η

0.5

η

Fig. 9. Dimensionless radial displacement (n = 1/2). Fig. 8. Dimensionless radial stresses (n = 0).

sical plate theory is still valid for the plate under consideration. It is further confirmed in Fig. 10, which shows that the axial displacement almost keeps invariant along the thickness direction. Thus a single deflection may be employed to represent the deformation of the plate, as adopted in the classical plate theory. Moreover, we can see that the deflection decreases with the gradient index k, since the plate becomes stiffer. The distributions of electric potential U and magnetic potential W along the thickness always look like a parabola, as shown in Figs. 11 and 12 respectively. It is also interesting to note that the magnetic potential has a butterfly-like shape of distribution that is very similar to that of the shear stress as shown in Fig. 3. When the plate is clamped, the distributions of rn, U, W, U and W have the same characters as those for the simply supported edge as shown above. We do not present them in the paper to save the space.

8 7 6

W (×10-2)

we see that the location of the maximum shear stress (here we means the absolute value jsngj) changes with the gradient index k. The point, at which the maximum shear stress is reached, approaches to the bottom surface (z = h/2) as k increases and the maximum shear stress increases with jkj. Fig. 4 shows that the distribution of electric displacement Dg in the homogeneous plate looks like a sinusoid. However, with the increase of jkj, it deviates from the sinusoidal shape and especially, the distribution curve for a positive k is quiet different from that for the negative k of the same absolute value. The same observation can be obtained when considering the distribution of Bg as shown in Fig. 6. The gradient index even has a more complicated influence on the distributions of Dn and Bn, as indicated in Figs. 5 and 7. Extensive numerical results indicate that rh(r,h/2) P rr(r,h/2) > 0 and rh(r,  h/2) 6 rr(r,  h/2) < 0 (0 6 r 6 a) always hold and the maximum hoop stress occurs at r = 0. The distribution of rh along the axisymmetric axis is the same as that of rr(rh = rr) and is depicted in Fig. 8 for the simply-supported conditions. It is seen that: (1) The maximum tensile and compressive stresses occur at (0, h/2) and (0, h/2), respectively; (2) The maximum tensile (compressive) stress increases (decreases) with the gradient index k. For example, the maximum tensile stress for k = 2.5 is less than the half of that in the homogeneous plate, while the maximum compressive stress is greatly reduced if k takes a value of 2.5. The distribution of the dimensionless radial displacement U at n = 1/2 with simply-supported edge is shown in Fig. 9, which indicates that the radial displacement changes always almost linearly along the thickness of the plate. This is somewhat different from the observation for electric displacement or magnetic induction as shown in Figs. 4–7, where the distribution shape greatly depends on the gradient index k. This implies that the assumption of the distributions of mechanical displacement in the clas-

k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

5 4 3 2 1 0 -0.5 -0.4

-0.3

-0.2

-0.1

0

η

0.1

0.2

0.3

Fig. 10. Dimensionless axial displacement (n = 0).

0.4

0.5

X.Y. Li et al. / Composite Structures 83 (2008) 381–390 4

2

Φ (×10-5)

0

-2 k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

-4

-6

-8 -0.5

-0.4

-0.3

-0.2

-0.1

0 η

0.1

0.2

0.3

0.4

0.5

Fig. 11. Dimensionless electric potential (n = 0). 4 2 k=-2.5 k=-1.0 k=0 k=1.0 k=2.5

0

-6 Ψ (×10 )

-2

389

(3) It is noteworthy that the direct displacement method can be applied to other axisymmetric problems. For example, if r5u5(z), r6w6(z), r4/4(z) and r4w4(z) are added to the expressions for u(r,z), w(r,z), /(r,z) and w(r,z) in Eqs. (7)–(9), respectively, analytical solution for a circular plate subjected to a load of parabolic form (qr2) can be obtained without extra difficulty. In fact, 3D solutions can be derived for loads with a general form of qr2n (n is a non-negative integer) simply by adding corresponding terms to the expressions for displacements, electric potential and magnetic potential in Eqs. (7)–(9). (4) Since the present solution was derived based on the three-dimensional theory of magneto–electro-elasticity without introducing any simplified assumptions on the elastic, electric and magnetic fields, it can serve as a benchmark for clarifying any approximate analysis or numerical method. Its importance is highlighted if we notice the fact that no numerical method performs well while keeping an economic computation cost, especially when the material properties of the plate vary rapidly through the thickness.

Acknowledgements

-4

The work was supported by the National Natural Science Foundation of China (Nos. 10472102 and 10432030). Partial support from the Specialized Research Fund for the Doctoral Program of Higher Education (No. 20060335107) is also acknowledged.

-6 -8 -10 -12 -0.5 -0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

η Fig. 12. Dimensionless magnetic potential (n = 0).

7. Conclusions (1) Three-dimensional analytical solution was derived for the axisymmetric problem of a uniformly loaded, transversely isotropic, circular plate of functionally graded magneto–electro-elastic materials. The material constants can vary along the thickness in an arbitrary way and the plate can be either simply-supported or clamped. The solution, which is independent of material eigenvalues [22], can be degenerated into that for homogeneous materials. (2) Numerical results show that the material inhomogeneity has an important effect on the elastic and electric fields in the plate. Thus, the plate can be designed by changing the material properties to meet any special requirements. In practice, however, the optimization object should be pre-determined case by case to make easy the choice of the spatial distribution of material. The explicitly expressed solution is very helpful to carry out such an optimization.

Appendix. Expressions for the elastic, electric and magnetic fields rz ¼ 64a1 H 1 ðzÞ þ 16a2 H 0 ðzÞ  q;

ðA1Þ

srz ¼ 8r½4a1 H 1 ðzÞ  a2 H 0 ðzÞ;

ðA2Þ

rr ¼ r2 ð4a1 z  a2 Þ½ð3c11 þ c12 Þ  4c13 f2 ðzÞ  4e31 g2 ðzÞ  4d 31 h2 ðzÞ þ 32a1 ½ðc11 þ c12 ÞF 1 ðzÞ  2c13 P 1 ðzÞ  2e31 Q1 ðzÞ  2d 31 R1 ðzÞ  8a2 ½ðc11 þ c12 ÞF 0 ðzÞ  2c13 P 2 ðzÞ  2e31 Q2 ðzÞ  2d 31 R2 ðzÞ  2a3 ½ðc11 þ c12 Þz  2c13 P 3 ðzÞ  2e31 Q3 ðzÞ  2d 31 R3 ðzÞ  2a4 ½ðc11 þ c12 ÞK 0 ðzÞ  2c13 P 4 ðzÞ  2e31 Q4 ðzÞ  2d 31 R4 ðzÞ  2a5 ½ðc11 þ c12 ÞK 1 ðzÞ  2c13 P 5 ðzÞ  2e31 Q5 ðzÞ  2d 31 R5 ðzÞ þ a6 ½ðc11 þ c12 Þ  2c13 P 6 ðzÞ  2e31 Q6 ðzÞ  2d 31 R6 ðzÞ  q½c13 P 7 ðzÞ þ e31 Q7 ðzÞ þ d 31 R7 ðzÞ;

ðA3Þ

390

X.Y. Li et al. / Composite Structures 83 (2008) 381–390

rh ¼ r2 ð4a1 z  a2 Þ½ðc11 þ 3c12 Þ  4c13 f2 ðzÞ  4e31 g2 ðzÞ  4d 31 h2 ðzÞ þ 32a1 ½ðc11 þ c12 ÞF 1 ðzÞ  2c13 P 1 ðzÞ  2e31 Q1 ðzÞ  2d 31 R1 ðzÞ  8a2 ½ðc11 þ c12 ÞF 0 ðzÞ  2c13 P 2 ðzÞ  2e31 Q2 ðzÞ  2d 31 R2 ðzÞ  2a3 ½ðc11 þ c12 Þz  2c13 P 3 ðzÞ  2e31 Q3 ðzÞ  2d 31 R3 ðzÞ  2a4 ½ðc11 þ c12 ÞK 0 ðzÞ  2c13 P 4 ðzÞ  2e31 Q4 ðzÞ  2d 31 R4 ðzÞ  2a5 ½ðc11 þ c12 ÞK 1 ðzÞ  2c13 P 5 ðzÞ  2e31 Q5 ðzÞ  2d 31 R5 ðzÞ þ a6 ½ðc11 þ c12 Þ  2c13 P 6 ðzÞ  2e31 Q6 ðzÞ  2d 31 R6 ðzÞ  q½c13 P 7 ðzÞ þ e31 Q7 ðzÞ þ d 31 R7 ðzÞ: Dr ¼ 2r½16a1 G1;z ðzÞ  4a2 G0;z ðzÞ  a4 T 11;z ðzÞ  a5 T 12;z ðzÞ; Dz ¼ 64a1 G1 ðzÞ þ 16a2 G0 ðzÞ þ 4a4 T 11 ðzÞ þ 4a5 T 12 ðzÞ: Br ¼ 2r½16a1 E1;z ðzÞ  4a2 E0;z ðzÞ  a4 T 12;z ðzÞ  a5 T 22;z ðzÞ; Bz ¼ 64a1 E1 ðzÞ þ 16a2 E0 ðzÞ þ 4a4 T 12 ðzÞ þ 4a5 T 22 ðzÞ; u ¼ rf4a1 ½8F 1 ðzÞ  r2 z  a2 ½8F 0 ðzÞ  r2   2a3 z  2a4 K 0 ðzÞ  2a5 K 1 ðzÞ þ a6 g; w¼

a1 ½r þ 16r f21 ðzÞ  64P 01 ðzÞ  4a2 ½r2 f20 ðzÞ  4P 02 ðzÞ þ a3 ½r2 þ 4P 03 ðzÞ þ 4a4 P 04 ðzÞ þ 4a5 P 05 ðzÞ  2a6 P 06 ðzÞ  qP 07 ðzÞ þ a8 :

/¼

16a1 ½r g21 ðzÞ  4Q01 ðzÞ  4a2 ½r2 g20 ðzÞ  4Q02 ðzÞ þ 4a3 Q03 ðzÞ þ 4a4 Q04 ðzÞ þ 4a5 Q05 ðzÞ  2a6 Q06 ðzÞ  qQ07 ðzÞ þ a9 ; 16a1 ½r2 h21 ðzÞ  4R01 ðzÞ  4a2 ½r2 h20 ðzÞ

w¼

4

ðA4Þ ðA5Þ ðA6Þ ðA7Þ ðA8Þ ðA9Þ

2

ðA10Þ

2

ðA11Þ

 4R02 ðzÞ þ 4a3 R03 ðzÞ þ 4a4 R04 ðzÞ þ 4a5 R05 ðzÞ  2a6 R06 ðzÞ  qR07 ðzÞ þ a0 :

ðA12Þ

where the involved functions can be found in Section 4 in the text. References [1] Tauchert TR, Ashida F, Noda N, Adali S, Verijenko V. Developments in thermopiezoelasticity with relevance to smart composite structures. Compos Struct 2000;48:31–8.

[2] Ootao Y, Tanigawa Y. Transient analysis of multilayered magneto– electro–thermoelastic strip due to nonuniform heat supply. Compos Struct 2005;68:471–80. [3] Li L, Zhou ZG, Wu LZ. Scattering of the harmonic elastic anti-plane waves from crack in piezoelectric/piezomagnetic materials. Chinese J Appl Mech 2006;23:101–4 [in Chinese]. [4] Zhou ZG, Chen Y, Wang B. The behavior of two parallel interface cracks in magneto–electro-elastic materials under an anti-plane shear stress loading. Compos Struct 2007;77:97–103. [5] Guan Q, He SR. Three-dimensional analysis of piezoelectric/piezomagnetic elastic media. Compos Struct 2007;77:97–103. [6] Harshe G, Dougherty JP, Newnham RE. Theoretical modeling of multiplayer magnetoelectric composites. Int J Appl Electromagnet Mater 1993;4:145–59. [7] Nan CW. Magnetoelectric effect in composites of piezoelectric and piezomagnetic phases. Phys Rev B 1994;50:6082–8. [8] Benveniste Y. Magnetoelectric effect in fibrous composites with piezoelectric and piezomagnetic phases. Phys Rev B 1995;51:16424–7. [9] Pan E. Exact solution for simply supported and multilayered magneto–electro-elastic plates. J Appl Mech 2001;68:608–18. [10] Zhou ZG, Wang B, Sun YG. Two collinear interface cracks in magneto–electro-elastic components. Int J Eng Sci 2004;42:1157–67. [11] Wu TL, Huang JH. Closed-form solutions for the magnetoelectric coupling coefficients in fibrous composites with piezoelectric and piezomagnetic phases. Int J Solids Struct 2000;37:2981–3009. [12] Wang X, Zhong Z. A finitely long circular cylindrical shell of piezoelectric/piezomagnetic composite under pressing and temperature change. Int J Eng Sci 2003;4:2429–45. [13] Mian AM, Spencer AJM. Exact solutions for functionally graded and laminated elastic materials. J Mech Phys Solids 1998;42(12):2283–95. [14] Li JY. Magneto–electro-elastic multi-inclusion and inhomogeneity problems and their applications in composite materials. Int J Eng Sci 2000;38:1993–2011. [15] Chen WQ, Lee KY. Alternative state space formulations for magnetoelectric thermoelasticity with transverse isotropy and the application to bending analysis of nonhomogeneous plate. Int J Solids Struct 2003;33:977–90. [16] Chen WQ, Lee KY, Ding HJ. On free vibration of non-homogeneous transversely isotropic magneto–electro-elastic plates. J Sound Vib 2005;279:237–51. [17] Pan E, Han F. Exact solution for functionally graded and layered magneto–electro-elastic plates. Int J Eng Sci 2005;43:321–39. [18] Bhangale RK, Ganesan N. Static analysis of simply supported functionally graded and layered magneto–electro-elastic plates. Int J Solids Struct 2006;43:3230–53. [19] Bhangale RK, Ganesan N. Free vibration of simply supported functionally graded and layered magneto–electro-elastic plates by finite element method. J Sound Vib 2006;294:1016–38. [20] Chen JY. Several problems of magneto–electro-elasticity. Doctorial dissertation, Zhejiang University, Hangzhou; 2004 [in Chinese]. [21] Ding HJ, Guo FL, Hou PF, Zou DQ. On the equilibrium of piezoelectric bodies of revolution. Int J Solids Struct 2000;37: 1293–326. [22] Ding HJ, Chen WQ. Three dimensional problems of piezoelasticity. New York: Nova Science Publishers; 2001.