Dynamic stability analysis of FGM plates by the moving least squares differential quadrature method

Composite Structures 77 (2007) 383–394 www.elsevier.com/locate/compstruct

Dynamic stability analysis of FGM plates by the moving least squares differential quadrature method Wu Lanhe *, Wang Hongjun, Wang Daobin Shijiazhuang Railway Institute, Department of Engineering Mechanics, Shijiazhuang 050043, China Available online 16 September 2005

Abstract The present paper investigates the dynamic stability of thick functionally graded material plates subjected to aero-thermomechanical loads, using a novel numerical solution technique, the moving least squares differential quadrature method. Temperature field is assumed to be a uniform distribution over the plate plane, and varied in the thickness direction only. Material properties are assumed to be temperature dependent and graded in the thickness direction in the simple power law manner. The equilibrium equations governing the dynamic stability of the plate are derived by the HamiltonÕs principle, then these equations are discretized by the moving least squares differential quadrature method. The boundaries of the instability region are obtained using the principle of BolotinÕs method and are conveniently represented in the non-dimensional excitation frequency to load amplitude plane. The influence of various factors such as gradient index, temperature, mechanical and aerodynamic loads, thickness and aspect ratios, as well as the boundary conditions on the dynamic instability region are carefully studied.  2005 Elsevier Ltd. All rights reserved. Keywords: Functionally graded plate; Dynamic stability; Forced vibration; Thermal load; Moving least squares differential quadrature method

1. Introduction The use of functionally graded materials has gained much popularity in recent years especially in extreme high temperature environments. Functionally graded materials (FGMs) are the new generation of composite materials in which the micro-structural details are spatially varied through non-uniform distribution of the reinforcement phase. This can be achieved by using reinforcement with different properties, sizes and shapes, as well as by interchanging the role of reinforcement and matrix phase in a continuous manner. The result is a macrostructure that produces continuous or smooth change on thermal and mechanical properties at the macroscopic level. Due to recent advances in material processing capabilities, that aid in manufacturing wide

*

Corresponding author. E-mail address: [email protected] (W. Lanhe).

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

variety of functionally graded materials, its use in application involving severe thermal environments is gaining acceptance in composite community, the aerospace and aircraft industry [1–4]. For instance, in a thermal protection system, FGMs take advantage of heat and corrosion resistance typical of ceramics and mechanical strength and toughness typical of metals. In view of these, there is an increased interest among researchers to study the dynamic behavior of the structural components made of these materials. It is easily found in the existing literature that many studies have been carried out on the vibration characteristics of isotropic plates and composite laminates. However, the investigations of buckling and vibration behaviors of functionally graded plates in thermomechanical environments are limited in number and addressed briefly herein. Tanigawa et al. [5] examined transient thermal stress distribution of FGM plates induced by unsteady heat conduction with temperature dependent properties. Reddy and Chin [6] dealt with

384

W. Lanhe et al. / Composite Structures 77 (2007) 383–394

many problems including transient response of FGM plate due to heat flux. Senthil and Batra [7] carried out 3-D analysis of transient thermal stresses in functionally graded plates adopting Laplace transformation technique and power series method. He et al. [8] presented a kind of finite element formulation based on thin plate theory for the shape and vibration control of FGM plate with integrated piezoelectric sensors and acuators under mechanical load whereas Liew et al. [9] have analyzed the active vibration control of plate subjected to thermal gradient using shear deformation theory. Yang and Shen [10] analyzed the dynamic response of initially stressed thin FGM plates subjected to impulsive load using GalerkinÕs procedure coupled with modal superposition method, and moreover, they studied the vibration characteristic and transient response of such plates in thermal environments by the first order shear deformation theory [11]. Cheng and Batra [12] studied the steady state vibration of a simply supported functionally graded polygonal plate with temperature independent material properties. Studies of static stability of FGM plates have received its due importance in the literature [13–16]. However, structural components, in general, under periodic loads can undergo parametric resonance that may occur over a range of forcing frequencies and if the load is compressive, resonance or instability can usually occur even if the magnitude of the load is below the static critical buckling load of the structure. So it is very important to study the dynamic stability behaviors of systems under periodic load [17–19]. Knowledge pertaining to such instability characteristics of functionally graded plate is meager in the literature and it is necessary for the optimal design and assessment of the structure failures. Ng et al. [20,21] presented the first contribution to the dynamic stability investigations of FGM structures based on the classical thin plate theory, without considering thermal environment and aerodynamic load. Furthermore, the governing equations were obtained using HamiltonÕs principle and the solutions were determined analytically employing the assumed mode technique. In practical applications, FGM plates are often used as structural components in flight structures, which are usually worked in the air and thermal environments. However, the influence of the airflow on the dynamic behaviors have not been studied in the exist reference, to the authorÕs knowledge. This is the motivation for this paper, the aim of which is to examine the effects of some factors such as temperature rise, gradient index, the airflow and etc., on the dynamic buckling behavior of FGM plates under periodic in-plane mechanical load. In the present paper, a novel numerical solution technique, the moving least squares differential quadrature method is employed to perform the study of the dynamic instability behavior of a thermo-mechanically

stressed FGM plate subjected to aerodynamic load. The temperature is assumed to be constant in the plane and varied only in the thickness direction of the plate. The material is assumed to be temperature dependent and graded in the thickness direction according to a power law distribution in terms of volume fractions of the constituents. A detailed investigation is carried out to bring out the influences of thickness and aspect ratios, thermal and mechanical loads, and boundary conditions on the dynamic instability region of functionally graded plates.

2. Material properties Consider a rectangular plate made of a mixture of metal and ceramic as shown in Fig. 1. The material in top surface and in bottom surface is metal and ceramic respectively. The modulus of elasticity E, the coefficient of thermal expansion a, the mass density q, the coefficient of thermal conductivity K and the PoissonÕs ratio m are assumed as Eðz; T Þ ¼ Ec ðT Þ  V c þ Em ðT Þ  ð1  V c Þ aðz; T Þ ¼ ac ðT Þ  V c þ am ðT Þ  ð1  V c Þ vðz; T Þ ¼ vc ðT Þ  V c þ vm ðT Þ  ð1  V c Þ

ð1Þ

qðzÞ ¼ qc V c þ qm ð1  V c Þ KðzÞ ¼ K c V c þ K m ð1  V c Þ where the subscripts m and c denote ceramic and metal respectively. Vc is the volume fraction of the ceramic and is assumed as a power function as follows V c ¼ ðz=h þ 1=2Þ

k

ð2Þ

where z is the thickness coordinate variable; and h/ 2 6 z 6 h/2, where h is the thickness of the plate and k is the power law index that takes values greater than or equal to zero. Here the coefficient of thermal conductivity K and the mass density q are assumed to be independent of temperature. Substituting Eq. (2) into Eq. (1), material properties of the FGM plate are determined, which are the same as the equations proposed by Praveen and Reddy [22]

a

b

y

x

h

z

Fig. 1. Configuration and coordinate system of a rectangular plate.

W. Lanhe et al. / Composite Structures 77 (2007) 383–394 k

where DT = Tc  Tm is defined as the temperature difference between ceramic-rich and metal-rich surface of the plate.

Eðz; T Þ ¼ Ecm ðT Þ  ðz=h þ 1=2Þ þ Em ðT Þ k

aðz; T Þ ¼ acm ðT Þ  ðz=h þ 1=2Þ þ am ðT Þ k

mðz; T Þ ¼ mcm ðT Þ  ðz=h þ 1=2Þ þ mm ðT Þ

385

ð3Þ

k

KðzÞ ¼ K cm ðz=h þ 1=2Þ þ K m

3. Basic equations

qðzÞ ¼ qcm ðz=h þ 1=2Þk þ qm where Ecm ðT Þ ¼ Ec ðT Þ  Em ðT Þ acm ðT Þ ¼ ac ðT Þ  am ðT Þ mcm ðT Þ ¼ mc ðT Þ  mm ðT Þ

ð4Þ

K cm ¼ K c  K m

Assume that u, v, w denote the displacements of the neutral plane of the plate in x, y, z directions respectively; hx, hy denote the rotations of the normal to the plate midplane. According to the first order shear deformation theory, the strains of the plate can be expressed as ex ¼ u;x þ zhx;x ; ey ¼ u;y þ zhy;y ;

qcm ¼ qc  qm

cxy ¼ u;y þ v;x þ zðhx;y þ hy;x Þ;

All the material properties of the two constituents such as E and a, which are temperature dependent, are assumed as

czx ¼ hx þ w;x ; czy ¼ hy þ w;y

P ¼ P 0 ðP 1 =T þ 1 þ P 1 T þ P 2 T 2 þ P 3 T 3 Þ

rx ¼

ð5Þ

where Pi(i = 1, 0, 1, 2, 3) are the coefficients of temperature T(K) and are unique for each constituent. For a functionally graded plate, the temperature change through thickness is not uniform, and is governed by the one-dimensional Fourier equation of heat conduction   d dT KðzÞ ¼ 0; T ¼ T c ðat z ¼ h=2Þ; dz dz T ¼ T m ðat z ¼ h=2Þ

"   kþ1 2z þ h K cm 2z þ h  2h 2h ðk þ 1ÞK m #  2kþ1 K 2cm 2z þ h þ 2h ð2k þ 1ÞK 2m "  3kþ1 DT K 3cm 2z þ h K 4cm  þ þ 3 C 2h ð3k þ 1ÞK m ð4k þ 1ÞK 4m  4kþ1  5kþ1 # 2z þ h K 5cm 2z þ h   2h 2h ð5k þ 1ÞK 5m

DT T ðzÞ ¼ T m þ C

ð7Þ with C ¼1 þ

K cm K 2cm K 3cm þ  ðk þ 1ÞK m ð2k þ 1ÞK 2m ð3k þ 1ÞK 3m K 4cm

ð4k þ 1ÞK 4m



HookeÕs law for a plate is defined as E ½ex þ mey  ð1 þ mÞaT ; 1  m2 E ½ey þ mex  ð1 þ mÞaT  ry ¼ 1  m2 E E E c ; szx ¼ c ; szy ¼ c sxy ¼ 2ð1 þ mÞ xy 2ð1 þ mÞ zx 2ð1 þ mÞ zy

ð10Þ The forces and moments per unit length of the plate expressed in terms of the stress components through the thickness are

ð6Þ

where Tc and Tm denote the temperature changes at the ceramic side and the metal side, respectively. The solution of Eq. (6) can be obtained by means of polynomial series [14]

K 5cm ð5k þ 1ÞK 5m

N ij ¼

Z

h=2

rij dz;

M ij ¼

h=2

Z

h=2

rij z dz; h=2

Qij ¼

Z

h=2

sij dz h=2

ð11Þ Substituting Eqs. (3), (9) and (10) into Eq. (11), gives the constitutive relations as E1 E2 ðu;x þ mv;y Þ þ ðhx;x þ mhy;y Þ  N Tx 2 1m 1  m2 E1 E2 Ny ¼ ðmu;x þ v;y Þ þ ðmhx;x þ hy;y Þ  N Ty 2 1m 1  m2 E1 E2 ðu;y þ v;x Þ þ ðhx;y þ hy;x Þ N xy ¼ 2ð1 þ mÞ 2ð1 þ mÞ Nx ¼

E2 E3 ðu;x þ mv;y Þ þ ðhx;x þ mhy;y Þ  M Tx 2 1m 1  m2 E2 E3 My ¼ ðmu;x þ v;y Þ þ ðmhx;x þ hy;y Þ  M Ty 2 1m 1  m2 E2 E3 ðu;y þ v;x Þ þ ðhx;y þ hy;x Þ M xy ¼ 2ð1 þ mÞ 2ð1 þ mÞ Mx ¼

Qx ¼ ð8Þ

ð9Þ

E1 ðhx þ w;x Þ; 2ð1 þ mÞ

Qy ¼

E1 ðhy þ w;y Þ 2ð1 þ mÞ ð12Þ

386

W. Lanhe et al. / Composite Structures 77 (2007) 383–394

where Z E1 ¼ E3 ¼

Z

h=2

Eðz; T Þ dz; h=2

h=2

Eðz; T Þz dz;

h=2

Eðz; T Þz2 dz

N Tx ¼ N Ty ¼ ¼

Z

h=2

h=2

M Tx

E2 ¼

E2 E2 ðu;xy þ v;xx Þ ðmu;xy þ v;yy Þ þ 2 1m 2ð1 þ mÞ E3 E3 ðhx;xy þ hy;xx Þ ðmhx;xy þ hy;yy Þ þ þ 1  m2 2ð1 þ mÞ E1 ðhy þ w;y Þ ¼ J h€y  2ð1 þ mÞ E1 E1 ðhx;x þ w;xx Þ þ ðhy;y þ w;yy Þ 2ð1 þ mÞ 2ð1 þ mÞ

M Ty

1 1m

1 ¼ 1m

Z

h=2

Eðz; T Þ  aðz; T Þ  T ðx; y; zÞ dz

h=2

Z

h=2

Eðz; T Þ  aðz; T Þ  T ðx; y; zÞz dz

h=2

qa V 2a ow ffi w € ¼q  N Tx w;xx  N Ty w;yy  N 0x w;xx þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M 21  1 ox

ð13Þ Suppose the plate is in the airflow environment, the aerodynamic pressure perpendicular to the plate surface q is determined by Birman and Librescu [23] qa V 2a ow ffi q ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M 21  1 ox

ð14Þ

If the plate is applied with an excitation load N 0x only in the x direction, the linear dynamic equations of equilibrium can be derived by the HamiltonÕs principle and are thus given by N x;x þ N xy;y ¼ 0 N y;y þ N xy;x ¼ 0 hx M x;x þ M xy;y  Qx ¼ J € €y M xy;x þ M y;y  Q ¼ J h

ð15Þ

y

w € Qx;x þ Qy;y  N Tx w;xx  N Ty w;yy  N 0x w;xx þ q ¼ q  and J are written as where q Z h=2 Z h=2 ¼ qðzÞ dz; J ¼ qðzÞz2 dz q h=2

ð17Þ By eliminating the variables u, v the equations of equilibrium can be covered into three equations in terms of hx, hy, w and are written as 1m E21 ð1  mÞ 1þm hx;yy  hy;xy hx þ 2 2 2 2ðE1 E3  E2 Þ E21 ð1  mÞ E1 ð1  m2 Þ € w  ¼ J hx ;x 2ðE1 E3  E22 Þ 2ðE1 E3  E22 Þ 1þm 1m E21 ð1  mÞ hx;xy þ hy;xx þ hy;yy  hy 2 2 2ðE1 E3  E22 Þ E21 ð1  mÞ E1 ð1  m2 Þ € w  ¼ J hy ;y 2ðE1 E3  E22 Þ 2ðE1 E3  E22 Þ 2ð1 þ mÞ T hx;x þ hy;y þ r2 w  ðN x w;xx þ N Ty w;yy þ N 0x w;xx Þ E1 2ð1 þ mÞ w € ¼ q ð18Þ E1 hx;xx þ

The general boundary conditions can be written as ð16Þ

h=2

When Eq. (12) are substituted into Eq. (14) the equations of equilibrium can be expressed in terms of displacement components. If the temperature distributes uniformly in x and y directions, the equations of equilibrium can be written as follows E1 E1 ðu;yy þ v;xy Þ ðu;xx þ mv;xy Þ þ 1  m2 2ð1 þ mÞ E2 E2 þ ðhx;yy þ hy;xy Þ ¼ 0 ðhx;xx þ mhy;xy Þ þ 2 1m 2ð1 þ mÞ E1 E1 ðu;xy þ v;xx Þ ðmu;xy þ v;yy Þ þ 1  m2 2ð1 þ mÞ E2 E2 ðhx;xy þ hy;xx Þ ¼ 0 ðmhx;xy þ hy;yy Þ þ þ 1  m2 2ð1 þ mÞ E2 E2 ðu;yy þ v;xy Þ ðu;xx þ mv;xy Þ þ 2 1m 2ð1 þ mÞ E3 E3 ðhx;yy þ hy;xy Þ ðhx;xx þ mhy;xy Þ þ þ 1  m2 2ð1 þ mÞ E1  ðhx þ w;x Þ ¼ J € hx 2ð1 þ mÞ

(1) Generalized hard simply supported edge (S) W ¼ 0;

M x ¼ 0;

hy ¼ 0

ð19Þ

(2) Clamped edge (C) W ¼ 0;

hx ¼ 0;

hy ¼ 0

ð20Þ

4. The MLSDQ method Consider a domain in the space X discretized by a set of discrete points {xi}i=1, 2, . . . , N. In the moving least squares differential quadrature (MLSDQ), the value of a partial derivative of a certain function u(x) at a discrete point xi can be approximate as a weighted linear sum of discrete function values chosen within the overall domain of a problem. hfuðxÞgi ¼

N X j¼1

cj ðxi Þuðxj Þ ¼

N X

cij uj

ð21Þ

j¼1

where denotes a linear differential operator which can be any orders of partial derivatives, cij are the weighting

W. Lanhe et al. / Composite Structures 77 (2007) 383–394

coefficients, and uj are the nodal function values. According to Liew et al. [24], the weighting coefficients cij can be determined directly from the partial derivatives of shape functions used in the moving least squares element-free method. Following Belytschko et al. [25], we have n X uh ðxÞ ¼ /i ðxÞui ð22Þ i¼1

where n is the number of nodes in the neighborhood of x and ui is the nodal parameter of u(x) at point xi, and ui the nodal shape function is given by 1

T

/i ðxÞ ¼ p ðxÞA ðxÞBi ðxÞ

pT ðxÞ ¼ ½1; x; y; x2 ; xy; y 2  n X  i ðxÞpðxi ÞpT ðxi Þ AðxÞ ¼ x B1 ðxÞ ½i¼1

BðxÞ ¼ B2 ðxÞ     i ðxÞpðxi Þ Bi ðxÞ ¼ x

ð24Þ ð25Þ Bn ðxÞ 

ð26Þ ð27Þ

 i ðxÞ ¼ xðx   xi Þ is a positive weighting function where x which decreases as jx  xij increases. It always assumes unity at the sampling point x and vanishes outside the domain of influence of x. The size of the domain of influence, or support size, determines the number of discrete points n in the domain of influence. Using Eq. (22), the weighting coefficients defined in the DQ method can be computed directly. For instance, we have h

ouðxÞ ou ðxÞ  ¼ ox ox

n X j¼1

discreted by a finite number of nodes (xi, yi), each of which are associated with three nodal parameters ðhix ; hiy ; wi Þ. Then, a circle of influence can be formed for each discrete point and the moving least-squares approximation of the displacement components are achieved in the domain of influence as described in above section hhx ðx; yÞ ¼

o/j ðxÞ uj ¼ ox

n X

ðxÞ

cj ðxÞuj

ð28Þ

j¼1

ðxÞ

where cj ðxÞ is the weighting coefficients of the first order derivative of u(x) in the x direction at any point x with respect to node xj. If the node xi is chosen as the sample point, the derivative at this node becomes   n n n X X X o/j ðxÞ ouðxÞ ðxÞ ðxÞ  u ¼ c ðx Þu ¼ cij uj  j i j j  ox x¼xi ox j¼1 j¼1 j¼1 x¼xi

ð29Þ Therefore, for each node xi in the domain, we have the corresponding weighting coefficient of the first derivative, ðxÞ such as cij , which are non-zero only for the points xj within the support domain of xi. The weighting coefficients for the other derivatives can be obtained in similar fashion.

5. Discretization of basic equations Now, letÕs discrete the governing equations and the boundary condition equations. First, the plate is

n X

/i ðx; yÞhix ; hhy ðx; yÞ ¼

i¼1

wh ðx; yÞ ¼

n X

n X

/i ðx; yÞhiy ;

i¼1

/i ðx; yÞwi

ð30Þ

i¼1

ð23Þ

where pi(x) is a finite set of basis functions. In this work on 2-D problems, the intrinsic polynomial basis is quadratic, i.e.

387

where n is the number of nodes within the domain of influence. It should be noted that n is dependent on the support size and may be different when different nodes are considered. One should also note that the nodal parameters are not equal to the physical values at the corresponding node. This is due to the shape functions used in the moving least squares method do not satisfy the KroneckerÕs delta condition generally, i.e., /i(xj, yj) 5 dij. Therefore, we have hx ðxi ; y i Þ  hhx ðxi ; y i Þ ¼

n X

/j ðxi ; y i Þhjx ¼

j¼1

hy ðxi ; y i Þ  hhy ðxi ; y i Þ ¼

n X

n X j¼1

/ij hjx

j¼1

/j ðxi ; y i Þhjy ¼

j¼1

wðxi ; y i Þ  wh ðxi ; y i Þ ¼

n X

n X

/ij hjy

ð31Þ

j¼1

/j ðxi ; y i Þwj ¼

n X

/ij wj

j¼1

where /ij = /j(xi, yi) for simplicity. Using the above equations together with the moving least squares differential quadrature technique, the governing Eq. (18) can be discretized at each point and, at the sampling point (xi, yi) these equations are written as 2

Ij11

Xn 6 6 6 0 j¼1 4 0

0 Ij22 0

2 0 n 6 X T 60 þ Nx 4 j¼1

0 8 9 0> > > < > = ¼ 0 > > > : > ; 0

38 j 9 2 j € > hx > > > R11 Rj12 > > > > n 6 7< = X 7 6 j j 0 7 € 6 R21 Rj22 hy > þ 5> 4 > > j¼1 > > > > j; Ij33 : w Rj31 Rj32 € 38 9 2 0 0 > hjx > 0 0 > > n < = X 7 6 0 j 60 0 0 0 7 5> hy > þ N x 4 > > j¼1 : ; j j 0 @33 0 0 w 0

38 9 Rj13 > hjx > = < > 7> j 7 R23 7 hjy > 5> > ; : j> Rj33 w 38 j 9 h > > x> > 7< j = 0 7 h 5> y > > : j> ; j }33 w 0

ð32Þ

where the coefficients Ijpp ðp ¼ 1; 2; 3Þ, Rjpq ðp; q ¼ 1; 2; 3Þ, @j33 ; }j33 are as follows

388

W. Lanhe et al. / Composite Structures 77 (2007) 383–394

E1 ð1  m2 ÞJ E1 ð1  m2 ÞJ /ij ; Ij22 ¼ /ij ; 2 2ðE1 E3  E2 Þ 2ðE1 E3  E22 Þ 2ð1 þ mÞ q Ij33 ¼ /ij ; E1 1  m ðyyÞ E21 ð1  mÞ ðxxÞ cij þ /ij ; Rj11 ¼ cij  2 2ðE1 E3  E22 Þ

inplane excitation. Assume that the excitation load is in the following form

Ij11 ¼

N 0x ¼ ðn þ f cos xtÞN 0cr

where N 0cr and x are the static buckling load and the excitation frequency of the forced inplane load. The dynamic instability region can be determined by the method suggested by Bolotin [17]. The principal instability region (first instability region) is usually the most important in studies of structures, both because of its widths as well as due to structural damping which often neutralizes higher regions. To obtain the boundaries of the instability region, the components {d} are written in the Fourier series as   1 X ixt ixt fdg ¼ þ fbgi cos fagi sin ð38Þ 2 2 i¼1;3;5;...

1 þ m ðxyÞ E21 ð1  mÞ ðxÞ cij ; Rj13 ¼ cij ; 2 2ðE1 E3  E22 Þ 1 þ m ðxyÞ c ; Rj21 ¼  2 ij 1  m ðxxÞ E21 ð1  mÞ ðyyÞ cij þ /ij ; Rj22 ¼ cij  2 2ðE1 E3  E22 Þ Rj12 ¼ 

Rj23 ¼

E21 ð1  mÞ ðyÞ cij ; 2ðE1 E3  E22 Þ ðyÞ

ðxÞ

Rj31 ¼ cij ; ðxxÞ

Rj31 ¼ cij ;

ðyyÞ

Rj33 ¼ cij  cij ; n 2ð1 þ mÞ X ðxxÞ ðyyÞ ðcij þ cij Þ; @j33 ¼  E1 j¼1

}j33 ¼ 

n 2ð1 þ mÞ X ðxxÞ cij E1 j¼1

ð33Þ

Similarly, the boundary condition equations can also be discreted into linear algebraic equations. (1) Generalized hard simply supported edge (S) n n X X /ij wj ¼ 0; /ij hjy ¼ 0; j¼1

E3 1  m2

j¼1 n X

ðxÞ

ðyÞ

ðmcij hjx þ cij hjy Þ  M Tx ðxi ; y i Þ ¼ 0

ð34Þ

j¼1

(2) Clamped edge (C) n n X X /ij hjx ¼ 0; /ij hjy ¼ 0; j¼1

j¼1

n X

/ij wj ¼ 0

j¼1

ð35Þ where i is the index number of the points on boundary. Combining all the discretized equations together with the corresponding boundary condition equations and rewritten them in terms of matrix form, one has ½Mf€ dg þ ð½K þ ½Kt þ N 0x ½Km Þfdg ¼ 0

ð37Þ

ð36Þ

where [M] and [K] are mass matrix and stiffness matrix, respectively; [K]t and [K]m are geometric stiffness matrix due to thermal and mechanical load, respectively; {d} is the displacement vector.

6. Instability analysis Eq. (36) represents the dynamic instability problem of a FGM plate subjected to thermal load and a parametric

with period 2T. This expression is substituted in Eq. (36) and the coefficients of each sine and cosine terms are set to be zero, as well as the sum of the constant terms. For nontrivial solutions, the determinants of the coefficients of these groups of linear homogeneous equations are equal to zero. Solving them for a given value of n, the variation of x with respect to f can be worked out. It has been shown by Bolotin that solutions with period 2T can be determined as a first order approximation from the equation       ½K þ ½Kt þ n  1 f ½Km  1 x2 ½M ¼ 0 ð39Þ   2 4 The solutions of Eq. (39) yield (1) static buckling factor under thermal load when f = x = 0, (2) free vibrations under thermal load when n = 0, f = 0, (3) free vibrations under thermal and axial constant mechanical load when f = 0, (4) regions of unstable solutions with constant n and varying f.

7. Numerical results and discussion The ceramic material used in this study is silicon nitride and the metal material used is stainless steel. The coefficient of thermal conductivity of these materials is assumed as constant, and they are taken to be 9.19 W/ mK and 12.04 W/mK respectively, for silicon nitride and stainless steel. The elastic moduli and the thermal expansion coefficients, as well as the PoissonÕs ratio however, are temperature dependent. The temperature coefficients of these materials are tabulated in Table 1. For convenience of presentation of the present results and for convenience of comparison with other numerical

W. Lanhe et al. / Composite Structures 77 (2007) 383–394

389

Table 1 Temperature dependent coefficients for silicon nitride Si3N4 and stainless steel SUS304 [27] Materials

Properties

P1

P0

P1 9

P2 4

P3 7

Si3N4

E (Pa) a (1/K) m q (kg/m3)

0 0 0 0

348.43 · 10 5.8723 · 106 0.24 2370

3.070 · 10 9.095 · 104 0 0

2.160 · 10 0.0 0 0

SUS304

E (Pa) a (1/K) m q (kg/m3)

0 0 0 0

201.04 · 109 12.33 · 106 0.3262 8166

3.079 · 104 8.086 · 104 2.002 · 104 0

6.534 · 107 0.0 3.797 · 107 0

results, a non-dimensional frequency parameter is defined as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 2 q ð1  m2 Þ m ð40Þ X¼x h Em Firstly, the present method is validated by considering the linear free vibration and static buckling analyses of FGM plates and by comparing the present results with other existing results. Table 2 presents the free vibration frequency parameters of FGM plates together with the numerical results given by Xiaolin and Huishen [26]. It is seen that the present results are in close agreement with those of Xiaolin. Table 3 presents the critical buckling temperature differences of such FGM plates together with the analytical solutions given by Lanhe

8.946 · 1011 0.0 0 0 0.0 0.0 0 0

[14]. It is also found that the present results agree very well with those of Lanhe. Based on the numerical comparison studies, we can confirm that the present method can yield accurate solutions for natural frequencies and buckling critical temperatures of functionally graded plates. The following study is focused mainly on the determination of boundaries of the primary instability regions that occur in the vicinity of simple resonance of first order. The plot of primary instability region in terms of non-dimensional excitation frequency parameter X versus the dynamic load amplitude f is depicted. Figs. 2 and 3 plot the dynamic stability results for a simply supported FGM plate under uniform temperature rise and no aerodynamic loading condition. It is seen that with increasing

Table 2 Comparison of frequency parameters of simply supported square FGM plate (a/h = 8) Temperature

k

Present

Xiaolin and Huishen [26]

Present

Xiaolin and Huishen [26]

Present

Xiaolin and Huishen [26]

Tc = 400, Tm = 300

0.0 0.5 1.0 2.0 10.0

12.353 8.513 7.439 6.678 5.742

12.397 8.615 7.474 6.693 –

29.033 20.115 17.578 15.717 13.558

29.083 20.215 17.607 15.762 –

43.775 30.226 26.510 23.699 10.602

43.835 30.530 26.590 23.786 –

Tc = 600, Tm = 300

0.0 0.5 1.0 2.0 10.0

11.958 8.253 7.144 6.378 5.423

11.984 8.269 7.171 6.398 –

28.433 19.468 17.116 15.375 13.138

28.504 19.784 17.213 15.384 –

43.003 29.886 25.994 23.315 20.083

43.107 29.998 26.109 23.327 –

X1

X2

X3

Table 3 Comparison of buckling temperature difference for simply supported square FGM plate (a/h = 10, Tm = 5) k

Resources

a/b = 1

a/b = 2

a/b = 3

a/b = 4

a/b = 5

0

Present Lanhe [14]

3258.427 3256.310

7645.309 7640.640

13840.534 13835.530

20771.455 20760.850

27598.918 27586.740

1

Present Lanhe [14]

1978.403 1976.297

4694.475 4691.691

8625.334 8619.424

13148.006 13141.430

17749.787 17740.250

5

Present Lanhe [14]

1482.105 1481.297

3480.570 3478.338

6293.531 6288.939

9420.568 9415.583

12486.072 12481.590

390

W. Lanhe et al. / Composite Structures 77 (2007) 383–394 0.7

0.7

k=0.0 k=0.5 k=5.0

0.5

0.6

Dynamic load amplitude

Dynamic load amplitude

0.6

0.4 0.3 0.2

k=0.0 k=0.5 k=5.0

0.5 0.4 0.3 0.2

0.1 0.1 0.0 12

16

20

24

28

0.0 8

Excitation frequency

12

16

20

24

Excitation frequency

Fig. 2. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/h = 20, n = 0.0, Tc = 300, Tm = 300, k = 0).

Fig. 4. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/h = 20, n = 0.2, Tc = 400, Tm = 300, k = 0).

0.7

0.7

k=0.0 k=0.5 k=5.0

0.5

0.6

Dynamic load amplitude

Dynamic load amplitude

0.6

0.4 0.3 0.2 0.1

k=0.0 k=0.5 k=5.0

0.5 0.4 0.3

0.2 0.1

0.0 12

16

20

24

28

Excitation frequency

Fig. 3. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/h = 20, n = 0.2, Tc = 300, Tm = 300, k = 0).

the value of gradient index k, the origin of instability shifts to lower forcing frequency, and the width of the instability region decreases for the given dynamic load amplitude f. The instability zone in general increases with dynamic load amplitude. It is also found that for a no mechanically pre-stressed plate with n = 0.0, the origin of the instability is higher than that for a prestressed plate with n = 0.2. Furthermore, it is observed that the boundaries of the instability regions with respect to graded index overlap at higher dynamic load, leading to wide range of operating frequency under which the system becomes unstable. Figs. 4 and 5 illustrate the dynamic instability results for a simply supported FGM plate with same temperature change Tc = 400, Tm = 300 and same static pre-stress parameter n = 0.2, respectively, with the

0.0 12

16

20

24

28

Excitation frequency

Fig. 5. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/h = 20, n = 0.2, Tc = 400, Tm = 300, k = 200).

aerodynamic load k = 200 and without aerodynamic load k = 0. It is observed that the presence of aerodynamic load can reduce the instability zone and postpone the occurrence of resonance. The origin of the instability changes to higher forcing frequency when aerodynamic load is present. Also, the instability regions for different gradient index k overlap at higher load amplitude compared to those of without airflow case. Figs. 6 and 7 compare the dynamic instability results for a simply supported FGM plate with different aspect ratio. It is obvious that with increasing of the plate aspect ratio, the origin of the instability region is changed to higher excitation frequency. It is also found that the width of the instability zone increases with the increase of aspect ratio.

W. Lanhe et al. / Composite Structures 77 (2007) 383–394 0.7

0.7

k=0.5 k=2.0 k=5.0

0.5

0.6

Dynamic load amplitude

Dynamic load amplitude

0.6

0.4

0.3 0.2

0.1

0.5

0.4

0.3 0.2

a/h=10 a/h=20 a/h=30

0.1

0.0

0.0 16

18

20

22

24

26

12

Excitation frequency

14

16

18

20

Excitation frequency

Fig. 6. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/h = 20, a/b = 1, n = 0.2, Tc = 400, Tm = 300, k = 400).

Fig. 8. Effect of length to thickness ratio a/h on the instability region for a simply supported FGM square plate (k = 0.5, a/b = 1, n = 0.2, Tc = 400, Tm = 300, k = 400).

0.7

0.7

0.6

0.6

k=0.5 k=2.0 k=5.0

0.5

Dynamic load amplitude

Dynamic load amplitude

391

0.4 0.3 0.2

0.5

0.4

0.3

0.2

a/h=10 a/h=20 a/h=30

0.1 0.1

0.0 20

24

28

32

36

40

44

48

Excitation frequency

Fig. 7. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/h = 20, a/b = 2, n = 0.2, Tc = 400, Tm = 300, k = 400).

Figs. 8 and 9 highlight the influence of the plate span to thickness ratio a/h on the dynamic instability region of a simply supported FGM plate for different gradient index. From these two Figures, it is seen that the thicker the plate, the higher the origin of the instability region will be. For the same gradient index k, the width of the instability region is reduced with increasing the plate thickness. For the low value of k, overlap in the instability regions occurs at relatively low excitation load amplitude when increasing the thickness. Figs. 10 and 11 depict the effect of temperature change on the dynamic instability results. Static mechanical load and aero-dynamic flow are also applied. It is seen that with the increase of the surface temperature difference, the origin of instability is shifted to

0.0 8

10

12

14

16

Excitation frequency

Fig. 9. Effect of length to thickness ratio a/h on the instability region for a simply supported FGM square plate (k = 5.0, a/b = 1, n = 0.2, Tc = 400, Tm = 300, k = 400).

lower forcing frequency and the instability frequency band will be more. Moreover, the overlap in the instability regions occurs relatively at low dynamic load amplitude while increasing the temperature difference. Figs. 12 and 13 focus on the effects of aerodynamic flow parameter k on the dynamic instability results of a simply supported FGM plate. It is viewed that when increasing the airflow parameter k, the resonance frequency will be increased for a given gradient index k. The widths of the unstable regions in the high airflow loading condition will be smaller than those in the low airflow loading condition, for a certain given index k. Figs. 14 and 15 describe the influence of boundary conditions on the dynamic instability results for a

W. Lanhe et al. / Composite Structures 77 (2007) 383–394 0.7

0.7

0.6

0.6

Dynamic load amplitude

Dynamic load amplitude

392

0.5

0.4

0.3

0.2

k=0.5 k=2.0 k=5.0

0.1

0.5

0.4 0.3

0.2

k=0.5 k=2.0 k=5.0

0.1

0.0

0.0 10

12

14

16

18

10

20

12

14

0.7

0.7

0.6

0.6

0.5

0.4

0.3 0.2

k=0.5 k=2.0 k=5.0

0.1

0.0 10

12

14

18

20

Fig. 12. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/b = 1, n = 0.2, Tc = 500, Tm = 300, k = 200).

16

18

Excitation frequency

Dynamic load amplitude

Dynamic load amplitude

Fig. 10. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/b = 1, n = 0.2, Tc = 400, Tm = 300, k = 200).

8

16

Excitation frequency

Excitation frequency

0.5

0.4 0.3

0.2

k=0.5 k=2.0 k=5.0

0.1

0.0 12

14

16

18

20

22

24

Excitation frequency

Fig. 11. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/b = 1, n = 0.2, Tc = 500, Tm = 300, k = 200).

Fig. 13. Effect of gradient index k on the instability region for a simply supported FGM square plate (a/b = 1, n = 0.2, Tc = 500, Tm = 300, k = 300).

FGM plate with the same thermal and mechanical load, for different gradient index. It is found that the origin and the band of instability regions are higher clamped boundary conditions than those for the simply supported boundary conditions, for both the cases k = 0.5 and k = 5.0.

ment and applied also with airflow load. Detailed formulations and analyses are presented. Several numerical examples are presented to demonstrate the effects of various factors such as gradient index, temperature difference, airflow load, plate aspect ratio, relative thickness and boundary conditions of the plate, on the dynamic instability results. From these examples, the following concluding remarks can be made:

8. Conclusions In this paper, the dynamic instability of functionally graded plates subjected to periodic inplane load is carefully studied by the moving least squares differential quadrature method. The plate is in the thermal environ-

(1) The origins and the widths of the unstable regions will be decreased with increasing the gradient index. For different gradient index, the dynamic instability regions will overlap at high amplitude of load.

W. Lanhe et al. / Composite Structures 77 (2007) 383–394

(4) With increasing of the surface temperature difference, the origin of instability is shifted to lower forcing frequency and the instability frequency zone will be more as expected. Moreover, the overlap in the instability regions occurs relatively at low dynamic load amplitude while increasing the temperature difference.

0.7

Dynamic load amplitude

0.6

0.5

0.4

0.3

References

0.2

0.1

Clamped Simply supported

0.0 12

14

16

18

20

22

24

26

28

30

32

34

Excitation frequency Fig. 14. Effect of boundary conditions on the instability region for a square FGM plate (a/b = 1, n = 0.2, Tc = 500, Tm = 300, k = 200, k = 0.5).

0.7

0.6

Dynamic load amplitude

393

0.5

0.4

0.3

0.2

0.1

Clamped Simply supported

0.0 8

10

12

14

16

18

20

22

24

Excitation frequency Fig. 15. Effect of boundary conditions on the instability region for a square FGM plate (a/b = 1, n = 0.2, Tc = 500, Tm = 300, k = 200, k = 5.0).

(2) The aerodynamic pressure load enhances the stability strength of the FGM plate by reducing the band of instability region and rising the forcing frequency of resonance. The higher the airflow load parameter k, the higher the value of resonance frequency will be. (3) When increases the aspect ratio and the thickness, as well as the boundary support conditions of the plate, the origin of instability changes to higher frequency, and the overlap of the instability regions occurs at relatively low dynamic load amplitude.

[1] Koizumi M. FGM activities in Japan. Composites 1997;28(1–2): 1–4. [2] Koizumi M. The concept of FGM. Ceramic Trans, Functionally Gradient Mater 1993;34:3–10. [3] Suresh S, Mortensen A. Functionally graded metals and metalceramic composites, Part 2: thermomechanical behavior. Int Mater Rev 1997;42:85–116. [4] Pindera MJ, Arnold SM, Aboudi J. Use of composites in functionally graded materials. Compos Eng 1994;4:1–145. [5] Tangiawa Y, Akai T, Kawamura R. Transient heat conduction and thermal stress problems of a nonhomogeneous plate with temperature-dependent material properties. J Therm Stresses 1996;19:77–102. [6] Reddy JN, Chin CD. Thermomechanical analysis of functionally graded cylinders and plates. J Therm Stresses 1998;21:593– 629. [7] Senthil SV, Batra RC. Three-dimensional analysis of transient thermal stresses in functionally graded plates. Int J Solids Struct 2003;40:7181–96. [8] He XQ, Ng TY, Sivashankar S, Liew KM. Active control of FGM plates with integrated piezoelectric sensors and actuators. Int J Solids Struct 2001;38:1642–55. [9] Liew KM, He XQ, Ng TY, Sivashankar S. Active control of FGM plates subjected to temperature gradient: modeling via finite element method based on FSDT. Int J Numer Meth Eng 2001;52:1253–71. [10] Yang J, Shen HS. Dynamic response of initially stressed functionally graded rectangular thin plates. Compos Struct 2001;54:497–508. [11] Yang J, Shen HS. Vibration characteristic and transient response of shear deformable functionally graded plates in thermal environments. J Sound Vibr 2002;255:579–602. [12] Cheng ZQ, Batra RC. Three dimensional thermoelastic deformations of a functionally graded elliptic plate. Compos, Part B: Eng 2000;31:97–106. [13] Esther F, Jacob A. Buckling analysis of functionally graded plates subjected to thermal loading. Compos Struct 1991;38: 29–36. [14] Lanhe W. Thermal bucking of a simply supported moderately thick rectangular FGM plate. Compos Struct 2004;64:211–8. [15] Najafizadeh MM, Eslami MR. First order theory based thermo elastic stability of functionally graded materials circular plates. AIAA J 2002;40:1444–50. [16] Ma LS, Wang TJ. Nonlinear bending and post-buckling of a functionally graded circular plate under mechanical and thermal loadings. Int J Solids Struct 2004;40:3311–30. [17] Bolotin VV. The dynamic stability of elastic systems. Sanfrancisco: Holden-Day; 1964. [18] Evan RM. On the parametric response of structures. Appl Mech Rev 1965;18:699–702. [19] Patel BP, Ganapathi M, Prasad KR. Dynamic instability of layered anisotropic composite plates on elastic foundations. Eng Struct 1999;21:988–95.

394

W. Lanhe et al. / Composite Structures 77 (2007) 383–394

[20] Ng TY, Lam KY, Liew KM. Effect of FGM materials on parametric response of plate structures. Comput Meth Appl Mech Eng 2000;190:953–62. [21] Ng TY, Lam KY, Liew KM, Reddy JN. Dynamic stability analysis of functionally graded cylindrical shells under periodic axial loading. Int J Solids Struct 2001;38:1295–309. [22] Praveen GN, Reddy JN. Nonlinear transient thermoelastic analysis of functionally graded ceramic-metal plates. Int J Solids Struct 1998;35:4457–76. [23] Birman V, Librescu L. Supersonic flutter of shear deformation laminated flat panel. J Sound Vibr 1990;139:265–75.

[24] Liew KM, Huang YQ. Bending and buckling of thick symmetric rectangular laminates using the moving least squares differential quadrature method. Int J Mech Sci 2003;45:95–114. [25] Belytschko T, Lu YY, Gu L. Element free Galerkin method. Int J Numer Meth Eng 1994;37:229–56. [26] Xiaolin H, Huishen S. Nonlinear vibration and dynamic response of functionally graded plates in thermal environments. Int J Solids Struct 2004;41:2403–27. [27] Yang J, Liew KM, Kitipornchai M. Dynamic stability of laminated FGM plates based on higher shear deformation theory. Comput Mech 2004;33:305–15.