Stability of piezoelectric FGM rectangular plates subjected to non-uniformly distributed load, heat and voltage

Advances in Engineering Software 39 (2008) 121–131 www.elsevier.com/locate/advengsoft

Stability of piezoelectric FGM rectangular plates subjected to non-uniformly distributed load, heat and voltage X.L. Chen b

a,b,*

, Z.Y. Zhao

a,b

, K.M. Liew

c

a Protective Technology Research Center, Nanyang Technological University, Nanyang Avenue, Singapore 639798, Singapore School of Civil and Environmental Engineering, Nanyang Technological University, Nanyang Avenue, Singapore 639798, Singapore c Department of Building and Construction, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong

Received 24 August 2005; received in revised form 15 December 2006; accepted 21 December 2006 Available online 27 February 2007

Abstract This paper employs the element free Galerkin (EFG) method to analyze buckling of piezoelectric functionally graded material (FGM) rectangular plates subjected to non-uniformly distributed loads, heat and voltage. A two-step solution procedure is implemented. In the first step, based on a two-dimensional (2D) plane stress problem, pre-buckling stresses of plates subjected to non-uniformly distributed loads are calculated. In the second step, based on the Mindlin plate assumption, buckling load and temperature parameters of the plates are obtained using the EFG method with penalty coefficients. Numerical examples are presented to validate the proposed numerical method. Some useful results for the FGM plates are presented.  2007 Elsevier Ltd. All rights reserved. Keywords: Functionally graded plate; Element free Galerkin method; Buckling analysis; Non-uniformly distributed load; Heat; Voltage

1. Introduction Composite materials with advantages of high stiffnessto-weight and strength-to-weight ratios are often used in aerospace engineering. They have drawn many researchers’ attention, and a lot of related literatures [1,2] can be found. However, composite materials have common limitations due to the feature of their layer-to-layer structures: delaminating between layers may occur due to the weakness of interfaces, debonding between matrix and fiber may exist due to the mismatch of mechanical properties when extreme thermal loads are applied, and residual stresses may present due to the difference in coefficients of thermal expansion of the fiber and the matrix. With the advent of

* Corresponding author. Address: Protective Technology Research Center, Nanyang Technological University, Nanyang Avenue, Singapore 639798, Singapore. Tel.: +65 6790 6895; fax: +65 6791 0676. E-mail address: [email protected] (X.L. Chen).

0965-9978/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2006.12.004

functionally graded materials (FGMs), the disadvantage of delaminating in composite materials can be overcome. The FGMs are a mixture of two or more materials in such way that their material volume fractions vary continuously along certain dimension. Due to the continuous variation of material properties from one surface to the other, the FGM plate avoids the interface problem that exists in homogeneous composite structures. The FGM can be often seen as thermal barrier plate structure that is subjected to high temperature. Japanese scientists first developed the FGM as early as in the late 1980s. Two decades later, the FGM has been widely used in several branches [3–7]. Buckling problem of FGM plates is very important and deserves us to pay more attention on it. Feldman and Aboudi [8] studied the elastic buckling solution for the FGM plate subjected to uniaxial loading. Liew et al. [9] analyzed the postbuckling of piezoelectric FGM plates subjected to thermo-electro-mechanical loading. Pan [10] gave the exact solution of simply supported and multilayered plates

122

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

subjected to magneto-electro-mechanical loading. It can serve as benchmarks to check various thick plate theories and numerical methods used for the modeling of layered composite structures. The plates that are subjected to non-uniformly distributed loads are often encountered in many engineering practices. However, solving buckling problem of such plates is not easy due to the existence of non-uniform distribution of pre-buckling stresses. Finite element method (FEM) has become a mainstream numerical method in engineering due to its versatility and the many powerful commercial software packages available. However, the existence of elements in the FEM sometimes embarrasses users. To avoid the use of elements, many sundry mesh-free methods have been developed recently. The key feature of the mesh-free methods is that the displacement approximation is implemented based on a set of scattered nodes instead of elements in the influence domain, thus the mesh-free methods can avoid the disadvantages that arise in the FEM from the use of elements, e.g., heavily distorted elements in the FEM may occur. Two kinds of meshfree methods deserve to be mentioned here. The element free Galerkin (EFG) method that employs scattered nodes without mass instead of elements to approximate displacements. The other one is the smoothed particle hydrodynamics (SPH) [11] that employs scattered particles with mass to approximate displacements. So far, many kinds of meshfree methods have been proposed and have their successful applications. Among the mesh-free methods, the element free Galerkin (EFG) method based on the moving least squares (MLS) approximation technique has the following advantages: (i) the shape function and its derivatives are smooth enough through the choice of appropriate weight function; (ii) the accuracy of the EFG results can reach high enough as desired for engineering application; (iii) the EFG results are very stable. For the plates subjected to non-uniform loads, the distribution of stresses in the plate domain changes largely. The solution of pre-buckling stresses of such plates is not easy for buckling problem. Seldom related literatures can be found. Liew and Chen [12–14] did some work on the buckling analysis of such plates. They employed the mesh-free method based on radial basis function and a two-step method to obtain buckling loads for rectangular and circular plates subjected to concentrated and partial uniform loads. This paper develops the EFG method for buckling analysis of the piezoelectric FGM plates that are subjected to non-uniformly distributed loads, electricity and heat. A two-step solution procedure is taken to implement the numerical computation. A variational form of static system equation is established based on a two-dimensional (2D) elastic plane stress problem in the first step, and the variational equation is discretized through displacement approximation. The pre-buckling stresses of the plane elasticity problem are obtained in this step when the plate is subjected to non-uniformly distributed loads. A variational

form of the system equation is established based on the Mindlin plate theory in the second step. The variational form is discretized to an eigenvalue equation. The static buckling loads and temperatures of such plates are obtained in this step. Buckling of the piezoelectric FGM plates, with or without electricity and temperature changes, and the pre-buckling stresses resulted from non-uniform loads, is studied.

2. Moving least squares (MLS) technique 2.1. Displacement approximation In this paper, the moving least squares (MLS) technique is used for displacement approximation. The elements and their connectivity used in the FEM are avoided. Only an array of scattered nodes in the domain under consideration is required for such approximation. Thus the approximation technique is flexible. The number and location of nodes can be handled freely and manually as user’s desire. Here we briefly introduce the generation of displacement approximation based on the MLS technique. The MLS approximation uh(x) of a displacement function u(x) at a point x in the problem domain is defined by a polynomial function as uh ðxÞ ¼

m X

pj ðxÞaj ðxÞ ¼ pt ðxÞaðxÞ

ð1Þ

j¼1

where p(x) is a complete basis whose m components are monomials in terms of space co-ordinates. a(x) is the vector of coefficients that correspond to the monomials and are functions of x. For buckling problem of plate structures, the dominant displacements are deflections caused by bending of plate structures. Hence we adopt a relatively high-order quadratic polynomial basis [15] as pt ðxÞ ¼ f1; x; y; x2 ; xy; y 2 g

ð2Þ

The unknown coefficients a(x) are obtained by minimizing a weighted residual function: J¼

 n X

~ ðx  xi Þ½pt ðxi ÞaðxÞ  ui 2 w

ð3Þ

i¼1

where n is the number of nodes in the domain under con~ ðx  xi Þ is a weight function sideration for the point x, w of compact support, ui is a nodal parameter. The stationary of J in Eq. (3) with respect to a(x) leads to the following linear algebraic equations: e e AðxÞaðxÞ ¼ BðxÞu s

ð4Þ

where us is the vector of nodal parameters of the displacements for all nodes in the support domain and is given by us ¼ fu1 ; u2 ; . . . ; un g

t

ð5Þ e and the symmetrical matrix AðxÞ and the non-symmetrical e matrix BðxÞ are defined by the following forms:

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

e AðxÞ ¼

n  X

t

 ¼ ~ ðdÞ ~ ðx  xi Þ ¼ w w



123

1  6d2 þ 8d3  3d4

for d 6 1 for d > 1

~ ðx  xi Þpðxi Þp ðxi Þ w

ð6Þ

e ~ ðx  xn Þpðxn Þ BðxÞ ¼ ½~ wðx  x1 Þpðx1 Þ; . . . ; w

ð7Þ

ð13Þ

Substituting a that is obtained from Eq. (4) into Eq. (1) yields

where d is the relative distance from node xi to point x that is defined by

0

i¼1

uh ðxÞ ¼

 n X

e 1 ðxÞ BðxÞÞ e ðpt ðxÞ A i ui ¼

i¼1

 n X

/i ðxÞui

ð8Þ

i¼1

where the vector of the shape functions is defined by UðxÞ ¼ f/1 ; /2 ; . . . ; /n g

ð9Þ

The shape functions satisfy linear consistency requirements: n X

/i ðxÞ ¼ 1

ð10aÞ

/i ðxÞxi ¼ x

ð10bÞ

/i ðxÞy i ¼ y

ð10cÞ

i¼1 n X i¼1 n X i¼1

jx  xi j d ¼ R0

ð14Þ

The value R0 can be chosen to be different for different point, but it is required that the number n of nodes in the support domain at any point x satisfies n P m. Thus, a e non-singular matrix AðxÞ in Eq. (4) can be ensured, and the unknown coefficients a in Eq. (4) can be obtained. We define R0 = dmaxc, in which dmax is a scaling parameter chosen by the users. c is determined at a particular node by searching for enough neighbor nodes to guarantee the e requirement that the matrix AðxÞ in Eq. (4) is non-singular and invertible everywhere in the support domain for computation of shape functions. 3. Governing equations

The shape functions do not satisfy the Kronecker delta criterion, namely

3.1. Thermo-electro-mechanical FGM plate

/i ðxj Þ 6¼ dij

A plate, shown in Fig. 1, is considered to be made of functionally graded materials (FGMs) through the thickness. The orthogonal Cartesian coordinates are introduced. The origin of the three axes is set in the undeformed plane of the plate. The FGM plate is considered to be made of combined ceramic–metal materials and its mechanical property is assumed to be linear elastic. The FGM plate varies smoothly and continuously through the thickness from pure metal of the top surface to pure ceramic of the bottom surface. The volume fractions of the un-symmetrically distributed materials are assumed to be [9]  n h þ 2z V c ðzÞ ¼ ð15aÞ 2h  n h þ 2z V m ðzÞ ¼ 1  ð15bÞ 2h

ð11Þ

where /i(xj) denotes the shape function of node xi within the support domain for the approximated point xj. Thus, one has uh ðxi Þ 6¼ ui

ð12Þ

Therefore, the MLS technique is an approximation rather than an exact displacement interpolation. This MLS property makes the imposition of essential boundary conditions more complicated than the FEM. 2.2. Weight function ~ ðx  xi Þ controls the size of supThe weight function w port domain and is very important for the performance of the MLS technique. The weight function should satisfy the following properties:

y,v

~ ðx  xi Þ should be positive and satisfy that 1. w jx  xij 2 R0, where R0 is the radius of compact support ~ ðx  xi Þ > 0 inside the support domain domain, i.e., w ~ ðx  xi Þ ¼ 0 outside the support domain. and w ~ ðx  xi Þ should be monotonically decreasing as the 2. w distance between point x and node xi increases. The function satisfying the above properties could be used as weight function. The quartic spline weight function and its first and second derivatives are continuous over the entire support domain and vanish at its support limit. For our application, we adopt the well-documented quartic spline weight function:

ha 2

h

ϕy

ha 2

x,u Metal

Ceramic

o

ϕx

b

Piezoelectric actuator

a z,w

Fig. 1. A FGM plate and its geometrical notation.

124

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

where Vc(z) and Vm(z) are the ceramic volume fraction and the metal volume fraction, respectively, and n is the nonnegative volume fraction exponent. The effective Young’s modulus E, Poisson’s ratio m and thermal expansion coefficient a can be expressed as  n h þ 2z E ¼ Em þ ðEc  Em Þ ð16aÞ 2h  n h þ 2z m ¼ mm þ ðmc  mm Þ ð16bÞ 2h  n h þ 2z a ¼ am þ ðac  am Þ ð16cÞ 2h For a plate that is subjected to non-uniformly distributed in-plane edge loads, pre-buckling stresses are non-uniform and cannot be determined using analytical method. Hence we adopt a two-step solution procedure: first, pre-buckling stress distribution in a plate is determined; second, buckling load of the plate with the pre-buckling stress distribution is obtained. For calculation of buckling temperature of the plate without applied in-plane edge loads, the first step is skipped.

According to Eq. (8) derived from the MLS technique, two components of the displacement vector u and v are approximated as u¼

 n X

/i ui

ð20aÞ

/i vi

ð20bÞ

i¼1

v¼

n  X i¼1

Substituting the approximated displacements u and v of Eqs. (20a) and (20b) into the variational formulation of Eq. (17) gives the discrete system equation as KU ¼ f

ð21Þ

where U is the vector of displacements of all nodes N in the entire problem domain and is defined by U ¼ fu1 ; v1 ; u2 ; v2 ; . . . ; uN ; vN gt

and K is the stiffness matrix that is assembled using the nodal stiffness matrix and is defined by [16] ! Z Z h=2 Kij ¼ Bti De Bj dz dA ð23Þ A

3.2. Prebuckling stress solution For the plate that is subjected to non-uniformly distributed in-plane edge loads, pre-buckling stresses can be calculated by taking the plate as a 2D elastic plane stress problem. The essential boundary conditions do not need to consider for the calculation of the pre-buckling stresses because the plate is balanced by the external in-plane edge loads. The variational formulation of the system is expressed as Z Z t de r dV  dutt dC ¼ 0 ð17Þ V

Cr

where e is the vector of strains with three components, r is the vector of stresses with three components, u is the vector of displacements with two components, and t is the vector of in-plane edge loads with two components per unit length. The formulations of the strain vector e and the stress vector r are 9 8 0 9 8 e > = > = < xx > < u;x > v;y ð18Þ e ¼ e0yy ¼ > > ; > ; : 0 > : u;y þ v;x cxy and 8 0 9 2 r 1 > = < xx > E 6 0 r r¼ ¼ 4m yy 2 > ; 1m : 0 > 0 sxy

38 0 9 e > = < xx > 7 e0 1 0 5 yy ¼ De e > ; : 0 > 0 ð1  mÞ=2 cxy m

0

ð19Þ where De is the elastic rigidity matrix, E is the effective elastic rigidity, and m is the effective Poisson’s ratio.

ð22Þ

h=2

and f is the force vector that is assembled using the nodal force vector and is defined by Z    /i t x fi ¼ dC ð24Þ  / Cr i ty in which Bi is defined by 2 3 /i;x 0 6 7 Bi ¼ 4 0 /i;y 5 /i;y

ð25Þ

/i;x

The displacements of the plate based on the 2D plane stress problem can be obtained by solving Eq. (21). Then, the stresses of any point can be obtained as follows: 8 0 9 >  n = X < rxx > r0yy ¼ De Bi ui ð26Þ > ; : 0 > i¼1 sxy where   ui ui ¼ vi

ð27Þ

3.3. Buckling solution of the thermo-electro-mechanical FGM plate For buckling analysis of a plate, the variational formulation with penalty coefficients of the plate can be written as Z Z Z t l t l t dðe Þ r dV þ den sn dV þ d ð~u  uÞ bð~u  uÞdC ¼ 0 V

V

Cu

ð28Þ

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

where el is the vector of linear strains, rl is the vector of linear stresses, en is the vector of non-linear strains, sn is the vector of stresses corresponding to en, b is the 5 · 5-dimensional and diagonal matrix with diagonal elements of big number defined by the users, ~ u is the boundary displacements, and  u is the prescribed boundary displacements. Based on the Mindlin plate assumption [17], the displacements of the plate can be written as uðx; y; zÞ ¼ u0 ðx; yÞ þ zux ðx; yÞ

ð29aÞ

vðx; y; zÞ ¼ v0 ðx; yÞ þ zuy ðx; yÞ

ð29bÞ

wðx; y; zÞ ¼ w0 ðx; yÞ

ð29cÞ

where u0, v0, and w0 are the displacements and transverse deflection of a point on the plate mid-plane, and ux and uy denote the rotations about the y-axis and the x-axis, respectively. Suppose that both the FGM and the piezoelectric material are linear elastic throughout the deformation. The linear strains el and the linear stresses rl can be expressed as elxx ¼ u0;x þ zux;x

ð30aÞ

elyy

¼ v0;y þ zuy;y

ð30bÞ

clxy

¼ u0;y þ v0;x þ zux;y þ zuy;x

ð30cÞ

clxz clyz

¼ w;x þ ux

ð30dÞ

¼ w;y þ uy

ð30eÞ

The linear stress–strain relationship for the FGM plate with piezoelectric material, taking into account the piezoelectric and thermal effects, is given by [9]

8 l 9 2 3 r > > 1 m 0 0 0 > > xx > > > > > 6 7 > > rl > > > 6m 1 7 0 0 0 > < yy > = 6 7 E 6 7 slxy ¼ 0 0 6 0 0 ð1  mÞ=2 7 2 > > 6 7 1  m > > > 60 0 7 l > > > 0 kð1  mÞ=2 0 sxz > > 4 5 > > > > > l ; > : 0 0 0 0 kð1  mÞ=2 syz 08 l 9 8 9 1 2 3 > a> 0 0 e31 > exx > > > > > > B> C > > > > 78 9 > elyy > > > >a> > C 6 B> 0 e32 7> Ex > > > > > > > > C 6 0 B> 7< = B< l = < = C 6 7 C 6 e15 0 0 7 Ey B B> cxy >  > 0 >DT C  6 > 6 7> > > > B> C 7: ; > cl > > > >0> > C 6 B> > > 4 0 e24 0 5 Ez xz > > > > > @> A > > > > > : cl > ; :0; 0 0 0 yz

ð31aÞ

125

where Va is the applied voltage through the thickness of the actuators. The high-order terms that are associated with the in-plane displacements in the non-linear strains en are neglected in this paper, and en can be simplified as 9 8 0 9 81 exx > > ðw;x Þ2 > > 2 > > > > > > > > > > > > > 1 ðw Þ2 > > e0 > > = = > < 2 ;y > < yy > 0 ð33Þ en ¼ cxy ¼ > > > > w;x w;y > > > > 0 > > > > > > > > > 0 > > cxz > > > ; ; > : : 0 > cyz 0 The stresses sn that corresponds to the non-linear strains en can be written as 8 0 9 8 0 9 rxx > rpxx > > > > > > > 0 > > 0 > > > > > > > > > > r > > < rpyy > < yy > = = 0 sn ¼ sxy ¼ kp s0pxy > > > > > > > > > > > 0 > > s0xz > > > > > > > > > > > : : 0 ; ; syz 0 2 3 1 m 0 0 0 6 7 0 0 0 6m 1 7 7 E 6 6 7 0 0 ð1  mÞ=2 0 0  kT 7 1  m2 6 6 7 0 kð1  mÞ=2 0 40 0 5 0 0 0 0 kð1  mÞ=2 8 9 2 3 0 0 e a > > 31 > > > > > 6 78 9 > >a> 0 0 e Ex > > > 6 32 7> < = 6 7< = 6 7 0 7 Ey ð34aÞ  0 DT  6 e15 0 > > : > ; > > 6 7> > > 0 e 0 E 0 > > 4 5 24 z > > > ; : > 0 0 0 0 and they can be rewritten in matrix form as  0 ¼ kp r 0p  kT DT Dl a  eE sn ¼ r

ð34bÞ

The five independent variables of Eqs. (29a)–(29c), u0, v0, w0, ux, and uy, of the plate can be approximated using Eq. (8). The shape functions for the five independent variables can be different. In this paper, for the simplicity of calculation the same set of shape functions is chosen for approximation of the five variables, which is expressed as u0 ¼

 n X

/i u0i

ð35aÞ

/i v0i

ð35bÞ

i¼1

and it can be expressed as rl ¼ Dl ðel  DT  aÞ  eE

ð31bÞ

where k is the shear correction factor (taken as 5/6 in Mindlin’s plate theory), Dl is the elastic rigidity matrix of the plate, and e31, e32, e15, e24 are the piezoelectric stiffness. The electric field is mainly dominated by its transverse component in the plate type piezoelectric material. Thus we can assume that f Ex

Ey

t

Ez g ¼ f 0

0

V a =ha g

t

ð32Þ

v0 ¼

 n X i¼1

w0 ¼

n  X

/i w0i

ð35cÞ

/i uxi

ð35dÞ

/i uyi

ð35eÞ

i¼1

ux ¼

 n X i¼1

uy ¼

 n X i¼1

126

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

Substituting the approximated independent variables of Eqs. (35a)–(35e) into the variational formulation of Eq. (28) gives the discrete eigenvalue equation for the buckling analysis of the plate as ½ðKl þ KE þ Ku Þ þ kp Gp þ kT GT Q ¼ RT þ RE

ð36Þ

where kp and kT are the eigenvalues corresponding to the applied load and the temperature change respectively, Q is the eigenvectors, and Kl, KE, Ku, Gp, and GT are the matrixes relating to the strain energy, the electric voltage, the penalty constraint, the applied load and the temperature change respectively. Kl is the matrix that is assembled using its nodal matrix Klij with 5 · 5 dimensions, i.e. ! Z Z h=2 t Klij ¼ Bli Dl Blj dz dA ð37Þ h=2

A

in which 2 /i;x 6 0 6 6 Bli ¼ 6 6 /i;y 6 0 4 0

0

0

z/i;x

/i;y /i;x

0 0

0 z/i;y

0 0

/i;x /i;y

/i 0

0

ð3;3Þ GTij

¼

Z "

ð/i;x /j;x þ /i;y /j;y Þ

Z

h=2 h=2

A

# aE dz DT dA ð1  mÞ ð43Þ

RT is the vector that is assembled by its nodal vector RTi with 5 · 1 dimensions, i.e. Z 1 RTi ¼ Bt Dl aDT dV ð44Þ 2 V li RE is the vector that is assembled by its nodal vector REi with 5 · 1 dimensions, i.e. Z 1 REi ¼ Bt eEdV ð45Þ 2 V li Solving Eq. (36), we can obtain the buckling loads or the buckling temperatures of the plate that is subjected to the combined thermal-electro-mechanical loadings.

3

z/i;y 7 7 7 z/i;x 7 7 0 7 5 /i

Table 1 Material properties of the FGM and piezoelectric materials

ð38Þ

KE is the matrix that is assembled using its nodal matrix KEij with 5 · 5 dimensions. ð3;3Þ KEij has the null elements except for the element KEij in third row and third column, i.e. Z ð3;3Þ ð39Þ KEij ¼ ha ½e31 /i;x /j;x þ e32 /i;y /j;y Ez dA

Material properties

Zirconia

Aluminum

G-1195N

Elastic modulus E (GPa) Poisson’s ratio m Thermal expansion coefficient a · 105 (1/C) Piezoelectric stress constant e31 (N/mV) Piezoelectric stress constant e32 (N/mV)

151.0 0.3 1.0

70.0 0.3 2.3

63.0 0.3 12.0

–

–

22.8429

–

–

22.8429

A

Ku is the matrix that is assembled using its nodal diagonal matrix Kuij with 5 · 5 dimensions, i.e. Z uti b/j dC ð40Þ Kuij ¼

a

Cu

b

in which 2 /i 6 60 6 ui ¼ 6 60 6 40 0

0

0

0

/i 0

0 /i

0 0

0

0

/i

0

0

0

0

3

7 7 7 7 7 7 05 /i 0 0

ð3;3Þ

Z

P

P

ð41Þ

Gp is the matrix that is assembled using its nodal matrix Gpij with 5 · 5 dimensions. ð3;3Þ Gpij has null elements except for the element Gpij in third row and third column: Gpij ¼ h

x

o

½r0pxx /i;x /j;x þ r0pyy /i;y /j;y þ s0pxy ð/i;x /j;y þ /i;y /j;x ÞdA A

ð42Þ GT is the matrix that is assembled using its nodal matrix GTij with 5 · 5 dimensions. ð3;3Þ GTij has null elements except for the element GTij in third row and third column:

y

Fig. 2. A rectangular plate subjected to a pair of concentrated loads.

Table 2 Buckling load parameter k cp ¼ PDcr0b of the simply supported plates subjected to a pair of concentrated loads Nodes

11 · 11

12 · 12

13 · 13

14 · 14

Isotropica (n = 0) FGM (n = 1)

26.527 19.019

25.682 18.630

25.572 18.525

25.573 18.539

a

Result of Leissa and Ayoub [18] is kcp = 25.814.

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

4. Numerical examples Rectangular plates of the FGM with a mixture of zirconia and aluminum are used for examination of the present EFG method. Material properties of the FGM and piezoelectric materials are listed in Table 1. The edges of the Table 3 Buckling temperature parameter kct = DTctac · 103 of the simply supported plates subjected to uniform temperature change

127

plates are a = b = 10.0 (m), and the thickness is h = 0.1 (m). The plates are subjected to non-uniformly distributed in-plane edge loads, heat and electricity. Two kinds of plates: Isotropic plate (i.e., volume fraction exponent n = 0) and FGM plate with volume fraction exponent n = 1 are considered. The buckling loads and temperatures of the simply supported and clamped plates are calculated using the present method. The boundary conditions are

-4

Nodes

11 · 11

12 · 12

13 · 13

14 · 14

8 x 10

Isotropica (n = 0) FGM (n = 1)

0.1288 0.08127

0.1266 0.07981

0.1261 0.07935

0.1264 0.07952

6

Result of Shen [19] is kct = 0.1264, and result of Noor and Burton [20] is kct = 0.1265.

4

a

SSSS, n=0 SSSS, n=1 CCCC, n=0 CCCC, n=1

Table 4 Buckling load parameter k cp ¼ PDcr0b of the plates subjected to a pair of concentrated loads SSSS

500.0 200.0 0.0 200.0 500.0

0 -2

CCCC

-4

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

25.677 25.615 25.573 25.531 25.469

18.644 18.581 18.539 18.498 18.435

69.997 69.958 69.932 69.905 69.862

49.061 49.019 48.987 48.954 48.900

-6 -8 -500 -400 -300 -200 -100

SSSS

500.0 200.0 0.0 200.0 500.0

100

200

300

400

500

Fig. 4. A square plate subjected to uniform temperature change.

a

Table 5 Buckling temperature parameter kct = DTctac · 103 of the plates subjected to uniform temperature change Va

0

Va

P

P

CCCC

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

0.1269 0.1266 0.1264 0.1262 0.1259

0.07998 0.07970 0.07952 0.07934 0.07906

0.3447 0.3443 0.3441 0.3439 0.3436

0.2100 0.2097 0.2096 0.2094 0.2091

c 2

c 2

b

Va

ΔKct

2

x

o c 2

c 2

0.15

y 0.1

SSSS, n=0 SSSS, n=1 CCCC, n=0 CCCC, n=1

ΔKcp

0.05

Fig. 5. A square plate subjected to a pair of partial uniform loads.

Table 6 Buckling load parameter k cp ¼ PDcr0b of the simply supported plates subjected to several pairs of partial uniform loads

0

-0.05

Va

-0.1 -0.15 -500 -400 -300 -200 -100

0

100

200

300

400

500

Va Fig. 3. A square plate subjected to a pair of concentrated loads.

500.0 200.0 0.0 200.0 500.0

c/b = 0.25

c/b = 0.5

c/b = 0.75

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

26.871 26.805 26.761 26.717 26.651

19.482 19.415 19.371 19.327 19.261

30.055 29.981 29.931 29.881 29.807

21.795 21.720 21.670 21.620 21.546

34.830 34.743 34.685 34.627 34.540

25.244 25.157 25.099 25.041 24.954

128

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

assumed to be u = v = w = 0 for simply supported edges and u = v = w = ux = uy = 0 for clamped edges. For calcu-

Table 7 Buckling load parameter k cp ¼ PDcr0b of the clamped plates subjected to several pairs of partial uniform loads Va

500.0 200.0 0.0 200.0 500.0

c/b = 0.25

c/b = 0.5

c/b = 0.75

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

73.474 73.413 73.373 73.332 73.272

51.458 51.398 51.357 51.316 51.255

81.701 81.632 81.586 81.540 81.472

57.218 57.150 57.104 57.058 56.989

92.165 92.087 92.034 91.982 91.903

64.546 64.467 64.415 64.362 64.284

c 2

c 2

P

P 2

P 2

x

b

o

a y

Fig. 6. A square plate subjected to axial and shear loads.

lation of buckling load, the temperature of the plate keeps unchanged. For calculation of buckling temperature, no in-plane edge loads are involved. For all examples considered here, the size of quadrilateral influence domain dmax [15] is chosen as 3.9 times the average nodal distance. Six terms’ complete polynomial bases are used. Regularly distributed nodes of 14 · 14 in the entire plate domain for displacement approximation and Gauss points of 4 · 4 in each background cell for numerical integration are used. 4.1. A pair of concentrated loads A rectangular plate that is subjected to a pair of in-plane concentrated loads, shown in Fig. 2, is considered. The non-dimensional buckling load parameter is defined as k cp ¼ PDcr0b, where Pcr is the critical loads that are applied to the edges of the plate, and the elastic rigidity is Ec h3 D0 ¼ 12ð1m 2 Þ. The non-dimensional buckling temperature parameter is defined as kct = DTctac · 103, where DTct is the critical temperature change, and ac is the coefficient of thermal expansion. First, the convergences of the buckling loads and temperatures of the isotropic and FGM simply supported plates are numerically studied using the proposed EFG method. Several densities of distributed nodes are used. The results are listed in Tables 2 and 3. When EFG nodes are less than 11 · 11, the results of FGM plate are not accurate enough. However, when EFG nodes are more than 11 · 11, the results fall in a reasonable range although the results for isotropic and FGM plates are a little oscillating due to employment of penalty coefficients in the EFG method. The convergence can be acceptable. The present

Table 8 Buckling load parameter k cp ¼ PDcr0b of the simply supported plates subjected to axial and shear loads Va

500.0 200.0 0.0 200.0 500.0

c/a = 0

c/a = 0.25

c/a = 0.5

c/a = 0.75

c/a = 1.0

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

32.251 32.198 32.162 32.126 32.072

23.051 22.999 22.963 22.928 22.875

34.066 34.008 33.969 33.930 33.872

24.365 24.308 24.270 24.231 24.174

38.904 38.835 38.789 38.743 38.674

27.859 27.791 27.746 27.700 27.632

46.528 46.443 46.386 46.329 46.244

33.368 33.284 33.228 33.172 33.087

56.017 55.911 55.840 55.768 55.662

40.212 40.106 40.035 39.965 39.859

Table 9 Buckling load parameter k cp ¼ PDcr0b of the clamped plates subjected to axial and shear loads Va

500.0 200.0 0.0 200.0 500.0

c/a = 0

c/a = 0.25

c/a = 0.5

c/a = 0.75

c/a = 1.0

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

84.755 84.704 84.670 84.637 84.586

59.454 59.403 59.369 59.335 59.284

89.825 89.770 89.734 89.697 89.642

62.989 62.934 62.898 62.861 62.806

103.328 103.263 103.220 103.177 103.112

72.442 72.377 72.334 72.291 72.226

123.707 123.626 123.572 123.518 123.437

86.716 86.635 86.581 86.527 86.446

147.984 147.881 147.813 147.744 147.641

103.719 103.616 103.547 103.479 103.376

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

converged results for isotropic plate are close to the other available results [18,19]. c 2

c 2 P

b

o

P 2

d

x

a

d

129

Second, the buckling loads and temperatures of the plate with the change of voltage Va are calculated. Their results are listed in Tables 4 and 5. To show the effect of voltage Va on the buckling parameters, the differences DKcp = Kcp  Kcp0 and DKct = Kct  Kct0, where Kcp0 and Kct0 are the buckling load and temperature parameters without voltage, i.e., Va = 0, respectively, are calculated and shown in Figs. 3 and 4, where SSSS denotes the simply supported plate while CCCC denotes the clamped plate. The buckling parameter differences arisen due to voltage for the isotropic and FGM plates with simply supported and clamped edges decrease linearly with a very small rate as the voltages increase as expected. The buckling parameters for the isotropic plate is bigger than those for the FGM plate, and the buckling parameters for the clamped plate is bigger than those for the simply supported plate. 4.2. Several pairs of partial uniform loads

P 2

y

A square plate that is subjected to several pairs of partial uniform in-plane edge loads, shown in Fig. 5, is considered.

Fig. 7. A square plate subjected to opposite axial loads.

Table 10 Buckling load parameter k cp ¼ PDcr0b of the simply supported plates subjected to opposite axial loads Va

c/a = 0

c/a = 0.25

d/a = 0

d/a = 0.25

d/a = 0

d/a = 0.25

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

500.0 200.0 0.0 200.0 500.0

28.898 28.842 28.805 28.767 28.711

20.727 20.671 20.634 20.598 20.542

29.376 29.318 29.279 29.240 29.182

21.072 21.015 20.977 20.939 20.881

30.175 30.115 30.075 30.035 29.976

21.660 21.601 21.562 21.523 21.463

30.700 30.638 30.596 30.555 30.493

22.036 21.975 21.934 21.893 21.831

500.0 200.0 0.0 200.0 500.0

33.549 33.480 33.434 33.388 33.320

24.113 24.045 24.000 23.954 23.886

34.172 34.100 34.052 34.004 33.932

24.560 24.489 24.441 24.394 24.322

43.745 43.648 43.584 43.519 43.422

31.520 31.423 31.359 31.294 31.198

44.272 44.169 44.100 44.031 43.927

31.903 31.800 31.732 31.663 31.560

c/a = 0.5

c/a = 1.0

Table 11 Buckling load parameter k cp ¼ PDcr0b of the clamped plates subjected to opposite axial loads Va

c/a = 0

c/a = 0.25

d/a = 0

d/a = 0.25

d/a = 0

Isotropic (n = 0)

FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

500.0 200.0 0.0 200.0 500.0

74.238 74.185 74.150 74.114 74.061

52.046 51.993 51.957 51.922 51.868

75.709 75.655 75.619 75.584 75.530

53.076 53.022 52.986 52.950 52.897

77.625 77.569 77.531 77.493 77.437

500.0 200.0 0.0 200.0 500.0

86.147 86.082 86.039 85.995 85.930

60.367 60.302 60.258 60.215 60.150

88.124 88.057 88.013 87.969 87.903

61.748 61.682 61.638 61.593 61.527

107.807 107.719 107.660 107.601 107.513

c/a = 0.5

Isotropic (n = 0)

d/a = 0.25 FGM (n = 1)

Isotropic (n = 0)

FGM (n = 1)

54.407 54.351 54.313 54.275 54.219

79.249 79.192 79.154 79.116 79.058

55.544 55.487 55.449 55.410 55.353

75.514 75.425 75.366 75.308 75.219

109.791 109.698 109.637 109.576 109.483

76.893 76.800 76.739 76.678 76.585

c/a = 1.0

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131

The applied edge loads vary from partial uniform loads to uniform loads through the change of c/b. The buckling loads of the isotropic and FGM plates that are subjected to several kinds of voltages are calculated using the present EFG method, and the results are listed in Tables 6 and 7. The buckling load parameters increase as the loads change from partial uniform loads to uniform loads. Again, the buckling load parameters decrease with a very small amount as the voltages increase.

110

d/a=0, n=0 d/a=0, n=1 d/a=0.25, n=0 d/a=0.25, n=1

100

90

k cp

130

80

70

4.3. Axial and shear loads A square plate that is subjected to one axial load and two shear loads, shown in Fig. 6, is considered. The axial load varies from concentrated load to uniform load through the change of c/a. The buckling loads of the plate that is subjected to several kinds of voltages are calculated. The results are listed in Tables 8 and 9. The buckling load parameters increase as the axial load changes from concentrated load to uniform load. The buckling load parameters decrease slightly as the voltages increase.

60

50

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

c/a (Va=-200.0 v) Fig. 9. Buckling load parameter k cp ¼ PDcr0b of the clamped plates subjected to opposite axial loads.

load parameters increase as the axial load on top edge changes from concentrated load to uniform load.

4.4. Three opposite axial loads 5. Conclusions A square plate that is subjected to three opposite axial loads, shown in Fig. 7, is considered. The axial loads vary through the change of c/a and d/a. The plate is subjected to several kinds of voltages. The calculated buckling load parameters are listed in Tables 10 and 11. Again, the buckling load parameters decrease with a small amount as the voltages increase. As an example, the buckling load parameters of the simply supported and clamped plate with the variation of axial load on the top edge and the fixed voltage Va = 200.0 V are calculated, and the results are drawn in Figs. 8 and 9. It can be found that the buckling

Acknowledgements

45

The first author would like to thank Prof. Liew KM because much part of the work of this paper was fulfilled when the first author was a research associate advised by Prof. Liew. The authors would also like to thank three reviewers for their valuable comments.

d/a=0, n=0 d/a=0, n=1 d/a=0.25, n=0 d/a=0.25, n=1

40

k cp

35

References

30

25

20

The calculated buckling load and temperature parameters are a little oscillating due to the employment of penalty coefficients in the proposed EFG method, but the results fall in a reasonable range, thus, the convergence still can be acceptable. The buckling parameters decrease with a small amount as the voltage increase. The buckling parameters for isotropic plates are bigger than those for FGM plates. The proposed EFG method and two-step solution procedure are valid for analysis of the buckling load and temperature of such piezoelectric FGM plates.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

c/a (Va=-200.0 v) Fig. 8. Buckling load parameter k cp ¼ PDcr0b of the simply supported plates subjected to opposite axial loads.

[1] Chen XL, Liu GR, Lim SP. An element free Galerkin method for the free vibration analysis of composite laminates of complicated shape. Compos Struct 2003;59:279–89. [2] Liu GR, Chen XL, Reddy JN. Buckling of symmetrically laminated composite plates using the element-free Galerkin method. Int J Struct Stab Dy 2002;2:281–94. [3] Liew KM, Kitipornchai S, Zhang XZ, Lim CW. Analysis of the thermal stress behaviour of functionally graded hollow circular cylinders. Int J Solids Struct 2003;40:2355–80. [4] Praveen GN, Reddy JN. Nonlinear transient thermoelastic analysis of functionally graded ceramic–metal plates. Int J Solids Struct 1998;35:4457–76.

X.L. Chen et al. / Advances in Engineering Software 39 (2008) 121–131 [5] Reddy JN. Analysis of functionally graded plates. Int J Numer Meth Eng 2000;47:663–84. [6] Reddy JN, Chin CD. Thermomechanical analysis of functionally graded cylinders and plates. J Therm Stresses 1998;26(1):93–126. [7] Woo J, Meguid SA. Nonlinear analysis of functionally graded plates and shallow shells. Int J Solids Struct 2001;38:7409–21. [8] Feldman E, Aboudi J. Buckling analysis of functionally graded plates subjected to uniaxial loading. Compos Struct 1997;38(1–4):29–36. [9] Liew KM, Yang J, Kitipornchai S. Postbuckling of piezoelectric FGM plates subjected to thermo-electro-mechanical loading. Int J Solids Struct 2003;40:3869–92. [10] Pan E. Exact solution for simply supported and multilayered magneto-electro-elastic plates. J Appl Mech-T ASME 2001;68: 608–18. [11] Liu GR, Liu MB. Smoothed particle hydrodynamics-a meshfree particle method. World Scientific; 2003. [12] Liew KM, Chen XL. Buckling of rectangular Mindlin plates subjected to partial in-plane edge loads using the radial point interpolation method. Int J Solids Struct 2004;41:1677–95. [13] Liew KM, Chen XL. Mesh-free radial point interpolation method for the buckling analysis of Mindlin plates subjected to in-plane point loads. Int J Numer Meth Eng 2004;60:1861–77.

131

[14] Liew KM, Chen XL, Reddy JN. Mesh-free radial basis function method for buckling analysis of non-uniformly loaded arbitrarily shaped shear deformable plates. Comput Method Appl M 2004;193: 205–24. [15] Liu GR, Chen XL. A mesh-free method for static and free vibration analyses of thin plates of complicated shape. J Sound Vib 2001;241:839–55. [16] Chen XL, Liew KM. Buckling of rectangular functionally graded material plates subjected to nonlinearly distributed in-plane edge loads. Smart Mater Struct 2004;13:1–8. [17] Mindlin RD. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. J Appl Mech-T ASME 1951;18: 31–8. [18] Leissa AW, Ayoub EF. Vibration and buckling of a simply supported rectangular plate subjected to a pair of in-plane concentrated forces. J Sound Vib 1988;127(1):155–71. [19] Shen HS. Thermal postbuckling analysis of imperfect Reissner– Mindlin plates on softening nonlinear elastic foundations. J Eng Math 1998;33:259–70. [20] Noor AK, Burton WS. Three-dimensional solutions for the thermal buckling of multilayered anisotropic plates. J Eng Mech-ASCE 1992;118:683–701.