Postbuckling analysis of functionally graded plates subject to compressive and thermal loads

Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

Contents lists available at ScienceDirect

Computer Methods in Applied Mechanics and Engineering j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / c m a

Postbuckling analysis of functionally graded plates subject to compressive and thermal loads Y.Y. Lee a,⁎, X. Zhao a, J.N. Reddy b a b

Department of Building and Construction, City University of Hong Kong, Kowloon Tong, Kowloon, Hong Kong Department of Mechanical Engineering, Texas A&M University, College Station TX, 77843-3123, USA

a r t i c l e

i n f o

Article history: Received 17 September 2009 Received in revised form 14 January 2010 Accepted 20 January 2010 Available online 1 March 2010 Keywords: Functionally graded materials Thermal and postbuckling analysis Functionally graded plates Ritz method von Kármán plates

a b s t r a c t Postbuckling analysis of functionally graded ceramic–metal plates under edge compression and temperature field conditions is presented using the element-free kp-Ritz method. The first-order shear deformation plate theory is employed to account for the transverse shear strains, and the von Kármán-type nonlinear strain– displacement relationship is adopted. The effective material properties of the functionally graded plates are assumed to vary through their thickness direction according to the power-law distribution of the volume fractions of the constituents. The displacement fields are approximated in terms of a set of mesh-free kernel particle functions. Bending stiffness is estimated using a stabilised conforming nodal integration approach, and, to eliminate the membrane and shear locking effects for thin plates, the shear and membrane terms are evaluated using a direct nodal integration technique. The solutions are obtained using the arc–length iterative algorithm in combination with the modified Newton–Raphson method. The effects of the volume fraction exponent, boundary conditions and temperature distribution on postbuckling behaviour are examined. © 2010 Elsevier B.V. All rights reserved.

1. Introduction The application of functionally graded materials (FGMs) in various engineering fields has been on the increase due to their special material properties, which change gradually and smoothly through certain dimensions according to a predetermined formula. FGMs are heterogeneous and typically made of metal and ceramic. Their unique features offer great potential for use in practical engineering applications, and have inspired considerable researcher interest. Numerous studies, involving thermal stress, vibration, buckling, static and dynamic analyses, have been carried out to date. Praveen and Reddy [1], for example, carried out transient thermoelastic analysis on functionally graded ceramic–metal plates subject to thermal and mechanical loads using a plate finite element method that accounts for the transverse shear strains, rotary inertia and moderately large rotations. Noda [2] discussed the thermal stress problems of FGMs, including the thermal stresses in a FGM plate, the thermal stress intensity factor of a FGM plate with a crack and transient thermoelasticity. Yang and Shen [3] investigated the free vibration and parametric resonance of shear deformable functionally graded plates (FGPs) in a thermal environment, and He et al. [4] conducted shape and ⁎ Corresponding author. Department of Building and Construction, City University of Hong Kong, Tatchee Avenue, Hong Kong Special Administrative Region, China. Tel.: + 852 2788 9847; fax: + 852 2788 7612. E-mail address: [email protected] (Y.Y. Lee). 0045-7825/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2010.01.008

vibration control analysis on FGM plates with piezoelectric sensors and actuators based on a finite element formulation. Liew et al. [5] investigated the static and dynamic piezothermoelastic controls of FGM plates embedded with integrated piezoelectric sensors and actuators in a thermal gradient environment. Efraim and Eisenberger [6] performed the vibration analysis of thick annular FGM plates of variable thickness. Reddy [7] developed a theoretical formulation and finite element model for a FGM plate using the third-order shear deformation theory, whereas Vel and Batra [8] reported a three-dimensional exact solution to the vibrations of FG rectangular plates. In addition to the aforementioned studies, researchers have also conducted various buckling and postbuckling analyses. For example, Najafizadeh and Heydari [9] investigated the thermal buckling of functionally graded circular plates based on the higher-order shear deformation plate theory, whereas Liew et al. [10] turned their attention to the postbuckling of FGM plates subject to thermoelectro-mechanical loads. Shen [11] presented thermal postbuckling analysis of FGM plates with temperature-dependent properties, and Woo et al. [12] developed an analytical solution to the postbuckling behaviour of moderately thick FGM plates. Using a 3-D finite element method, Na and Kim [13–15] conducted three dimensional thermomechanical buckling and thermal postbuckling analysis of FGM plates. The research methods commonly adopted in the postbuckling analysis of structures include analytical [16–18], semi-analytical [19] and finite element methods [20–22]. As analytical methods are

1646

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

suitable only for simple problems, and the solutions obtained from finite element methods depend strongly on the quality of the meshes, mesh-free methods, in which the interpolations are based entirely on discrete nodes rather than meshes, have recently been developed and are now often adopted as an alternative approach in engineering analysis. The element-free Galerkin method, for example, has been employed to study static and dynamic fractures [23], thin plates [24] and thin shells [25]. The meshless local Petrov–Galerkin method has been applied to thin beam analysis [26], and the reproducing kernel particle method has been employed in numerical simulations of the large deformations of thin shell structures [27] and in metal-forming analysis [28]. The spline strip kernel particle method has been used to analyse folded plate structures [29]. In this paper, the postbuckling behaviour of FGM plates under edge compression and temperature field conditions are investigated using the element-free kp-Ritz method [30–32]. The first-order shear deformation plate theory and the Von Kármán assumption are employed to account for transverse shear strains, rotary inertia and moderate rotations, and the material properties are assumed to vary continuously and smoothly through the thickness according to the power-law distribution of the volume fractions of the constituents. The temperature field is considered to vary through the thickness direction alone and to be held constant in any plane. The bending stiffness of the plates is assessed using a stabilised conforming (SC) nodal integration approach [33], and the shear and membrane terms are obtained employing a direct nodal integration method. Compared with Gauss integration, the proposed integration scheme improves computational efficiency and avoids shear locking for very thin plates. To solve the nonlinear system equations, the arc–length method [34] is employed in combination with the modified Newton-Raphson technique to track the complete postbuckling path. The current formulation is verified by comparisons between the present results and those given in the literature, and the influence of the volume fraction exponent, side-to-thickness ratio, temperature field and boundary conditions on the postbuckling behaviour of FGM plates is discussed. 2. Functionally graded plates Fig. 1 shows a ceramic–metal FGP measuring a, b and h (length, width and thickness). A coordinate system (x, y, z) is established on the middle plane of this plate. The material properties are assumed to

Fig. 2. The volume fraction versus the thickness.

vary through the thickness according to the following power law distribution. P ðzÞ = ðPc −Pm ÞVc + Pm ;

Vc =


 1 z n + ðn≥0Þ; 2 h

ð1aÞ ð1bÞ

where P represents the effective material properties, including Young's modulus E, density ρ, Poisson's ratio v, thermal conductivity k and thermal expansion α; Pc and Pm denote the properties of the ceramic and the metal, respectively; Vc is the volume fraction of the ceramic; and n is the volume fraction exponent. Fig. 2 describes the variation in the volume fraction through the thickness for different volume fraction exponents n. For temperature-dependent materials, their properties are given by   −1 2 3 P = P0 P−1 T + 1 + P1 T + P2 T + P3 T ;

ð2Þ

where P0, P− 1,P1, P2 and P3 are the temperature coefficients. In this study, the material properties are computed at room temperature (T = 27 °C) unless otherwise specified. 3. Theoretical formulation 3.1. Energy functional According to the first-order shear deformation plate theory [35], the displacement field can be written as uðx; y; zÞ = u0 ðx; yÞ + zϕx ðx; yÞ vðx; y; zÞ = v0 ðx; yÞ + zϕy ðx; yÞ wðx; y; zÞ = w0 ðx; yÞ;

ð3Þ

where u0, v0 and w0 denote the displacements at the mid-plane of the plate in the x, y and z directions, and ϕx and ϕy represent the transverse normal rotations about the y and x axes, respectively. The nonlinear strain–displacement relationship is given by

Fig. 1. A functionally graded plate.

8 9   < εxx = γyz εyy = ε0 + zκ; = γ0 ; γxz :γ ; xy

ð4a; bÞ

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

K = 5/6 for isotropic materials; for functionally graded materials, it is modified as [6]

where 9 > > > > > > > > > = ; ε0 = > > > > > > > > > > > > > ∂u ∂v0 ∂w0 ∂w0 > > > > > + ; : 0 + ∂y ∂x ∂x ∂y 8 > > > > > > > > > <


 ∂u0 1 ∂w0 2 + 2 ∂x ∂x
 ∂v0 1 ∂w0 2 + 2 ∂y ∂y

K= ð5Þ

9 > > 9 8 > > > ∂w0 > > > > > > > > > + ϕ = < y ∂y = κ= ; γ0 = : > > > ∂y > > > > ∂w0 > > > > > > ; > : + ϕ > x > > ∂ϕy > > > ∂x ∂ϕ > > ; : x + ∂y ∂x 8 > > > > > > > > <

∂ϕx ∂x ∂ϕy

ð6a; bÞ

2

Q11 6 Q12 6 6 0 6 4 0 0

Q12 Q11 0 0 0

0 0 Q66 0 0

0 0 0 Q44 0

9 8 9 1 308 1> 0 εxx > > > > > > > > > > > > > > > B C 7 0 7B< εyy = < 1 = C B γxy − 0 αΔT C; ð7Þ 0 7 C 7B> > > > A 0 5@> γ > > >0> > > > > ; : > ; > : yz > γxz 0 Q55

where ΔT is the temperature change from a stress-free state (T = 0 °C), and the elastic coefficient Qij is a function of z because the effective properties of the plate vary along the thickness direction according to Eqs. (1a, 1b). The total in-plane force resultants, total moment resultants and transverse force resultants are given by 8 9 8 9 8 9 8 9 < Nxx = < σxx = < Mxx = < σxx = h= 2 h= 2 N = Nyy = ∫−h = 2 σyy dz; M = Myy = ∫−h = 2 σyy zdz; :N ; :σ ; :M ; :σ ; xy xy xy xy ð8a; bÞ  Qs =

Qy Qx





h=2

= ∫−h = 2

 σyz dz: σxz

ð9Þ

N M



   ε   ˆ  A B 0 − N ; =  κ ˆ B D M s

Q s = A γ0 ;

ð16Þ

1 ˜ T εdΩ; ∫ N 2 Ω

ð17Þ

where Ñ and ɛ are defined as 8 9 8 9 < ε0 = < N = ˜ = κ : N M ;ε= : ; : ; γ0 Qs

ð18a; bÞ

The external work carried out by the surface traction and the body force is written as T

We = ∫Ω u

 T  fdΩ + ∫Γ u tdΓ;

ð19Þ

_ _ where f and t represent the external load and prescribed traction, respectively, on the natural boundary. Thus, the total potential energy functional for the plate can be expressed as Πs = Uε −We :

ð20Þ

3.2. Two-dimensional kernel particle shape functions

ð10Þ

The construction of the kernel particle shape functions is briefly reviewed here. Consider a two-dimensional domain, Ω, that is discretized by a set of scattered nodes xI, I = 1, ⋯⋯, NP. The approximation of displacements is denoted by uh and expressed in terms of the kernel particle shape functions as

ð11Þ

u = ∑ ψI ðxÞuI ;

In matrix form, the relation between the stress resultants and the strains can be written as 

5 ; 6−ðv1 V1 + v2 V2 Þ

where V1 and V2 represent the volume fractions of two different components, and ν1 and ν2 are the Poisson's ratios of the two constituents. This modified shear correction coefficient accounts for the effects of different constituents in terms of their volume fractions and Poisson's ratios, therefore yields more accurate results than the constant coefficient. It has been demonstrated that different shear correction coefficients have no pronounced influence on the solutions of functionally graded plates when the length to thickness ratio is greater than 10 [36]. The strain energy of the plate is given by Uε =

The constitutive relations are expressed as 9 8 σ > > > > > xx > > = < σyy > σxy = > > > σ > > > > ; : yz > σxz

1647

h

NP

ð21Þ

I=1

where N̂ and M̂ represent the thermal force and moment resultants, respectively, and are given by

where uI represents the nodal coefficient, and ψI(x) denotes the shape function associated with node I and can be defined as [37,38]

ˆ = ∫h = 2 ½ 1 N −h = 2

1

0  ðQ11 + Q12 ÞαΔTdz;

ˆ = ∫h = 2 ½ 1 M −h = 2

1

0  ðQ11 + Q12 ÞαΔTzdz:

T

ð12Þ

ψI ðxÞ = C ðxÞΦa ðx−xI Þ;

T

ð13Þ

where Φa(x − xI) is the kernel function, and C(x) is the correction function and is used to satisfy the reproducing condition:

_ The matrices A, B, D, and As are the inplane, bending–stretching coupling, bending, and shear stiffnesses [35], and their corresponding terms, Aij, Bij, Dij and Asij, are given by     h= 2 2 Aij ; Bij ; Dij = ∫−h = 2 Qij 1; z; z dz;

ð14Þ

NP

p q

p q

∑ ψI ðxÞxI yI = x y

for

ð22Þ

ð23Þ

p + q = 0; 1; 2:

I=1

Correction function C(x) is expressed as a combination of complete second-order monomial functions: T

C ðxÞ = H ðx−xI ÞbðxÞ; s Aij

=

h= 2 K∫−h = 2 Qij dz;

ð15Þ Asij

is defined for i, where Aij, Bij and Dij are defined for i, j = 1, 2, 6, and j = 4, 5. The transverse shear correction coefficient is taken to be

ð24aÞ T

bðxÞ = ½b0 ðx; yÞ; b1 ðx; yÞ; b2 ðx; yÞ; b3 ðx; yÞ; b4 ðx; yÞ; b5 ðx; yÞ ; T

2

2

H ðx−xI Þ = ½1; x−xI ; y−yI ; ðx−xI Þðy−yI Þ; ðx−xI Þ ; ðy−yI Þ ;

ð24bÞ ð24cÞ

1648

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

where H is a vector with a quadratic basis and b(x) is a coefficient vector that has yet to be determined. The shape function can now be written as T

ψI ðxÞ = b ðxÞHðx−xI ÞΦa ðx−xI Þ:

ð25Þ

The coefficient b(x) can be solved by substituting Eq. (25) into Eq. (23), as follows. −1

bðxÞ = M

ðxÞHð0Þ;

ð26Þ

where M is termed a moment matrix and H(0) is a constant vector. The expressions for M and H(0) are given by NP

T

MðxÞ = ∑ Hðx−xI ÞH ðx−xI ÞΦa ðx−xI Þ;

ð27aÞ

I=1

T

Hð0Þ = ½1; 0; 0; 0; 0; 0 :

ð27bÞ

Φa ðx−xI Þ = Φa ðxÞ⋅Φa ðyÞ;

ð29Þ

where φ(x) is the weight function. In this study, the cubic spline function is chosen as the weight function, and is given by

where zI =

ðx−xI Þ , dI

9 > > > > > = ; > > > > > ;

ð30Þ

dI is the size of the support and is given by

dI = dmax cI ;

where dmax is a scaling factor, and distance cI is chosen such that there is a sufficient number of nodes to avoid the singularity of matrix M. Finally, the shape function is expressed as T

−1

ðxÞHðx−xI ÞΦa ðx−xI Þ:

where ɛij is the strain and AL is the area of representative domain ΩL. Applying the divergence theorem to Eq. (33) yields the following strain smoothing expression. ε˜ ij ðxL Þ =

  1 h h ∫ u n + uj ni dΓ; 2AL ΓL i j

ð34Þ

where ΓL represents the boundary of the representative domain of node L and n is the outward normal of boundary ΓL. The discrete form of the strain smoothing is given by ð35Þ

2

˜ ðx Þ b I1 L 6 ˜ BI ðxL Þ = 4 0 ˜ ðx Þ b I2 L

3 0 ˜ ðx Þ 7 5; b I2 L ˜b ðx Þ I1 L

˜ ðx Þ = 1 ∫ ψ ðx Þn ðx ÞdΓ b Ii L AL ΓL I L i L

ð36Þ

ði = 1; 2Þ ;

ð37Þ

where SL is a group of nodes with a shape function that supports cover node L. Eq. (37) can be estimated using Gauss integration or the trapezoidal rule. Note that no derivatives of the shape functions are required for the integration: the two-dimensional domain integration is simply transferred to the one-dimensional line integration, which significantly reduces the computational cost. 4. Discrete system equations

ð31Þ

ψI ðxÞ = H ð0ÞM

ð33Þ

I∈SL

in which

1 for 0≤jzI j≤ 2 1 bjz j≤1 for 2 I otherwise

1 ∫ ε ðxÞdΩ; AL ΩL ij

h

ð28Þ

x−x  I ; Φa ðxÞ = φ a

ε˜ ij ðxL Þ =

˜ ðx Þu ; ˜ε ðxL Þ = ∑ B I L I

The kernel function Φa(x − xI) is expressed as

8 2 > 2 3 > −4zI + 4zI > > > <3 4 4 3 φz ðzI Þ = 2 > −4zI + 4zI − zI > > 3 3 > > : 0

Consider a domain, Ω, that is discretised by a set of discrete nodes, xI, I = 1, ⋯⋯, NP, and the representative domain of node xL is denoted by ΩL, which can be generated from the Voronoi diagram shown in Fig. 3. The strain smoothing at node xL is defined as

4.1. Thermal analysis For a plate subject to a temperature field, which is assumed to vary only along the thickness direction, the temperature distribution along

ð32Þ

As the shape function, ψI(x), does not possess Kronecker delta properties, the boundary conditions cannot be directly imposed. Several methods, such as the transformation method [37], Lagrange multipliers [24] and the penalty method, are available to enforce the essential boundary conditions. 3.3. Nodal integration There are two integration schemes available for mesh-free methods to evaluate the stiffness matrix: Gauss quadrature and nodal integration. The former is computationally cost-ineffective, whereas the latter is more efficient but not very stable due to the reduced integration [39,40]. To address these shortcomings, Chen et al. [33] proposed SC nodal integration, which ensures linear exactness and eliminates the spurious oscillations inherent in direct nodal integration. A number of researchers have demonstrated the efficiency and accuracy of this integration method in dealing with solid mechanics problems [41–43]. The essence of SC nodal integration is the utilisation of strain smoothing over the representative integration domain. A brief review is provided here.

Fig. 3. The nodal representative domain obtained from a Voronoi diagram.

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

1649

Matrices KL, KNL, u and F are given by L

b

m

sh

t

K =K +K +K

+K;

ð43Þ

u = ½ u1

T

u2

⋯ un  ;

ð44Þ

b

bT

b

ð45Þ

m

mT

m

sh

sT

s s

KIJ = ∫Ω BI DBJ dΩ; mT

KIJ = ∫Ω BI ABJ dΩ + ∫Ω BI

 b bT  m BBJ dΩ + ∫Ω BI BBJ dΩ;

ð46Þ

KIJ = ∫ Ω BI A BJ dΩ; T ˆ  t KIJ = ∫Ω GI N GJ dΩ NL

KIJ = ∫Ω


 1 LT NL 1 NLT NL NLT L BI SBJ + BI SBJ + BI SBJ dΩ 2 2

" " # ˆ   N m T b T FI = ∫Ω ψI fdΩ + ∫Γ ψI tdΓ + ∫Ω BI BI  dΩ; ˆ M

Fig. 4. Temperature profiles through the plate thickness.

that direction is governed by the one-dimensional steady-state heat conduction equation, as follows.
 d dT kðzÞ = 0: − dz dz

ð47a; bÞ

ð38Þ

2

m

BI

ð48Þ

ð49Þ

3

6 b 7 NL  L 7 BI = 6 4 BI 5; BI = H G;

ð50a; bÞ

BsI 2

The boundary conditions are T = Ttop

at z = h = 2 ;

ð39aÞ

T = Tbot

at z = −h = 2 ;

ð39bÞ

where Ttop and Tbot represent the temperatures on the top and bottom surfaces, respectively, and k(z) denotes the thermal conductivity. Fig. 4 describes the temperature variation with various volume fraction exponents for an Al/ZrO2 plate; the top and bottom surfaces are held at temperatures of Ttop = 200 °C and Tbot = 20 °C, respectively. Using the temperature distribution results, the thermal resultants, N̂ and M̂ , can be calculated according to Eqs. (12) and (13).

3  A B 0  4 S = B D 0 5: 0 0 As

ð51Þ

The stiffness matrix in Eq. (45) is evaluated using the stabilised nodal integration technique, and Eqs. (46)–(49) are computed by direct nodal integration. The approximations of Eqs. (45)–(49) are given as follows. NP

bT

˜ ðx ÞDB ˜ ðx ÞA ; KIJ = ∑ B J I L L L b

b

ð52Þ

L=1

4.2. Discrete system equations

 NP  b m mT m mT KIJ = ∑ BI ðxL ÞABJ ðxL Þ + BI ðxL Þ BBJ ðxL Þ L=1  m bT + BI ðxL Þ BBJ ðxL Þ AL ;

The discrete form of the approximations of the displacements of the middle plane of a plate can be expressed as

KIJ = ∑

0

1 uh0 0 1 B C uI B vh C B 0C B vI C B C NP NP B C h B C C u0 = B wh0 C = ∑ ψI B ψI ðxÞuI : B wI C = I∑ B C I=1 =1 @ ϕxI A B hC B ϕx C ϕ @ A

ð40Þ

NL

NP



1 LT L NLT L B ðxL ÞSBJ ðxL Þ + BI ðxL ÞSBJ ðxL Þ 2 I 1 NLT NL + BI ðxL ÞSBJ ðxL Þ AL ; 2 L=1

NP NP T  sh sT s s t KIJ = ∑ BI ðxL ÞA BJ ðxL ÞAL ; KIJ = ∑ GI ðxL Þ N GJ ðxL ÞAL ; L=1

ð53Þ

L=1

ð54Þ

ð55Þ

yI

h

ϕy

Substituting Eq. (40) into Eq. (20) yields the following discrete system equation. Ks ðuÞu = F;

ð41Þ

where stiffness matrix Ks is expressed as L

NL

Ks ðuÞ = K + K ðuÞ;

ð42Þ

in which KL represents the linear stiffness matrix and KNL denotes the nonlinear stiffness matrix, which is a function of the displacements.

NP NPb  FI = ∑ ψI ðxL ÞfðxL ÞAL + ∑ ψI ðxL Þ tðxL ÞsL L=1 L =1   NP ˆ mT bT + ∑ BI ðxL Þ BI ðxL Þ N AL ; ˆ L=1 M

ð56Þ

where xL and AL are the node coordinate and nodal representative area, respectively, NPb is the number of nodes on the natural boundary, and sL are the weights associated with the boundary s ̅ ̅ point. Here, B̃bI (xL), BbI (xL), Bm I (xL), BI (xL), H and G are given by 2

3 ˜ ðx Þ 0 0 0 0 b Ix L ˜Bb ðx Þ = 6 ˜ ðx Þ 7 60 0 0 7 0 b I L Iy L 5; 4 ˜ ˜ 0 0 0 bIy ðxL Þ bIx ðxL Þ

ð57Þ

1650

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

2 60 0 0 6 6 6 b BI ðxL Þ = 6 60 0 0 6 6 4 0 0 0

3

∂ψI ðxL Þ ∂x

0

7 7 7 ∂ψI ðxL Þ 7 7; ∂y 7 7 7 ∂ψI ðxL Þ 5

0 ∂ψI ðxL Þ ∂y

ð58Þ

∂x

˜ ðx Þ = 1 ∫ ψ ðxÞn ðxÞdΓ; ˜ ðx Þ = 1 ∫ ψ ðxÞn ðxÞdΓ; b b Ix L x Iy L y AL ΓL I AL ΓL I 2

∂ψI ðxL Þ 6 ∂x 6 6 6 m BI ðxL Þ = 6 6 0 6 6 4 ∂ψI ðxL Þ ∂y 2 60 0 s 6 BI ðxL Þ = 6 4 0 0

ð59a; bÞ

3 07 7 7 7 0 0 07 7; 7 7 5 0 0 0

0

0

∂ψI ðxL Þ ∂y ∂ψI ðxL Þ ∂x

∂ψI ðxL Þ ∂y

0

∂ψI ðxL Þ ∂x

ψI ðxL Þ

0

ð60Þ

5.1. Postbuckling analysis of isotropic plates

3 ψI ðxL Þ 7 7 7; 5 0

ð61Þ

2

3T ∂w ∂w 0 0 07 6 ∂x ∂y  6 7 H=6 7 ; 4 5 ∂w ∂w 0 0 0 ∂y ∂x 2 60 0  6 G=6 4 0 0

∂ψI ðxL Þ 0 ∂x ∂ψI ðxL Þ 0 ∂y

Fig. 5. The geometry of a plate subject to edge compression.

To verify the present formulation, an isotropic square plate subject to edge compression, as shown in Fig. 5, is considered first. Yamaki [17] investigated this case using an analytical method, and Dawe and Wang [44] employed a spline finite strip method to do so. The lengthto-thickness ratio is a/h = 120, and the Poisson's ratio is v = 1/3. All of the edges are simply supported, and the unloaded edges are free to move laterally along the plane of the plate. The boundary condition is expressed as

ð62Þ SSSS-01: at loaded edges x = 0, a at unloaded edges y = ±b/2

3  07 N 7  7; N = 0xx 5 0

0 : Nyy

ð63a; bÞ

In this study, the arc–length iterative strategy is combined with the modified Newton–Raphson method [34], in which the tangent stiffness matrix is computed only at the beginning of each load step and then held constant throughout the subsequent iterative cycles, and employed to determine the entire load–deflection path. 5. Numerical examples and discussion In this section, the postbuckling responses of FGPs under conditions of edge compression and temperature field are investigated. The FGPs considered here comprise aluminium and zirconia. The properties of each constituent, including the Young's modulus, Poisson's ratio, thermal expansion coefficient, conductivity and density, are given in Table 1. The kernel particle function is employed to construct the shape functions of the displacement expressions, and a scaling factor of 2.2, which represents the size of the support, is adopted. The stiffness matrices are estimated using the nodal integration method described in Section 3.3, and the essential boundary conditions are enforced using the transformation method. To trace the complete postbuckling path, a small initial deflection in the transverse direction is imposed on the plates, and the modified Newton–Raphson method is combined with the arc–length approach to solve the nonlinear equilibrium equations.

v = w = ϕy = 0; w = ϕx = 0.

Non-dimensional load factor N̆ is defined as N̆ = Nxa2/π2Eh3, where E is the elastic modulus, and Nx is the uniformly distributed load along the edges of x = 0, a. Fig. 6 shows the variation in non-dimensional load factor N̆ against non-dimensional central displacement w/h. The present results are obtained using a nodal distribution of 13 × 13, and the two sets of results, namely, Present-01 and Present-02, are produced using different initial deflections. For comparison, the solutions given in References [17,44] are also plotted in this figure. Note that the results labelled Yamaki-01 and Yamaki-02 were obtained for plates without and with initial deflections, respectively. It can be seen that the results achieved using the proposed methods agree well with those reported in the literature. Next, a square plate with the same geometrical and material properties as those in the last case, but with a different boundary condition, is considered. The

Table 1 Properties of the FGM components. Material

Aluminium (Al) Zirconia (ZrO2)

Properties E (N/m2)

v

α (/°C)

κ (W/Mk)

70.0 × 109 151 × 109

0.30 0.30

23.0 × 10− 6 10.0 × 10− 6

204 2.09

Fig. 6. The postbuckling response of a simply supported square plate (SSSS-01) subject to edge compression.

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

1651

Fig. 7. The postbuckling response of a square plate (CCSS) subject to edge compression.

Fig. 9. The postbuckling response of aluminium/zirconia FGM plates subject to edge compression, with all edges simply supported (SSSS-02, a/h = 100).

loaded edges of this plate are clamped, and the unloaded edges are simply supported. The boundary condition is expressed as follows.

N̑ = Nxa2/Emh3, where Em represents the Young's modulus of the metal, and Nx is the uniformly distributed load along the edges of x = 0, a. All of the results in this section are obtained using 15 × 15 nodes, and the volume fraction exponents are taken as n = 0, 0.5, 1, 2, ∞. Fig. 8 depicts the postbuckling responses of simply supported (SSSS-02) aluminium/zirconia plates with a thickness ratio of a/h = 50 and subject to edge compression. It can be seen that, in general, the pure zirconia and aluminium plates achieve the highest and lowest load factors, respectively. The load factors of the plates with n = 0.5, n = 1 and n = 2 fall between those of the pure zirconia and aluminium plates. When n is closer to 0, the load factor of the plate is closer to that of the pure zirconia plate. This is as expected because the stiffness of the plates decreases as the volume fraction increases. A similar phenomenon can also be discerned in Figs. 9 and 10, which describe the postbuckling behaviour of aluminium/zirconia plates with a/h = 100 and a/h = 200, respectively. Fig. 11 illustrates the postbuckling responses of aluminium/zirconia plates with the CCSS boundary condition and a thickness ratio of a/ h= 100. A similar trend to that seen in Fig. 8 is observed, except that the corresponding load factor is higher than that of the plate with the SSSS-02 boundary condition. This is because more constraints are imposed on the plates with the CCSS boundary condition. Fig. 12 depicts the influences of the two simply supported boundary conditions (SSSS-01 and SSSS-02) on the postbuckling responses of the aluminium/zirconia plates. It is

CCSS: at loaded edges x = 0, a at unloaded edges y = ±b/2

v = w = ϕx = ϕy = 0; w = ϕx = 0.

The present results are compared with those given by Yamaki [17], as shown in Fig. 7, in which the curves labelled Yamaki-01 and Yamali-02 were obtained for a plate with and without initial deflections, respectively. Again, good agreement can be observed. 5.2. Postbuckling analysis of FGM plates In this subsection, the postbuckling analysis of aluminium/zirconia FGM plates subject to edge compression is examined. Two types of boundary conditions, named SSSS-02 and CCSS, are considered. The CCSS condition is as defined in the previous case, and SSSS-02 is expressed as follows. SSSS-02: at loaded edges x = 0, a at unloaded edges y = ±b/2

v = w = ϕy = 0; u = w = ϕx = 0.

The geometrical properties of the FGM plates used here are a/ b = 1, a/h = 50, 100, 200. Load factor parameter N̑ is defined as

Fig. 8. The postbuckling response of aluminium/zirconia FGM plates with all edges simply supported (SSSS-02, a/h = 50) under the condition of edge compression.

Fig. 10. The postbuckling response of aluminium/zirconia FGM plates with all edges simply-supported (SSSS-02, a/h = 200) and subject to edge compression.

1652

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

Fig. 11. The postbuckling response of aluminium/zirconia FGM plates (CCSS, a/h = 100) subject to edge compression.

Fig. 13. The postbuckling response of aluminium/zirconia FGM plates edge compression and temperature field conditions (SSSS-02, a/h = 20).

found that, at the initial stage when the deflection is small, the postbuckling curves of the two are only slightly different; when the load factor increases further, however, the plates with the SSSS-01 boundary condition undergo larger deformations than those with the SSSS-02 boundary condition. The effects of temperature on the postbuckling responses of FGM plates are also investigated. A simply supported (SSSS-02) FGM aluminium/zirconia plate subject to edge compression and a temperature field is considered. The geometrical properties of this plate are a =b = 0.2 m and h = 0.01 m, and the volume fraction exponent is taken as n = 0.5. The bottom surface of the plate is held at Tbot = 20 °C, and the top surface is held at Ttop = 100 °C; the temperature distribution can be solved according to Eq. (38). Fig. 13 shows the postbuckling behaviour of the plate subject to edge compression and a temperature field. For comparison, the postbuckling curve of the plate under the edge compression condition alone is also plotted. It can be seen that, due to the imposed temperature field, the plate experiences larger postbuckling deformations than in the case without the temperature field. This is because the temperature field produces a thermal force and moment resultants and reduces the stiffness of the plate. The effects of different temperature fields on the postbuckling responses of FGM plates are also examined. The aforementioned plate is used here, but three different temperature fields are imposed upon it: the bottom surface remains at Tbot = 20 °C, whereas the temper-

ature on the top surface varies from 100 °C to 300 °C. Fig. 14 reveals the postbuckling responses of the plate subject to these three temperature fields. A similar trend can be observed in the three postbuckling curves, and the plate subject to the highest temperature field undergoes the largest initial deflection and a greater degree of postbuckling deformation, which indicates that a higher temperature field leads to a reduction in the stability of the plate and a greater degree of postbuckling deformation.

The postbuckling responses of FGM plates under edge compression and temperature field conditions have been investigated using the kp-Ritz method. The effective material properties of FGPs are assumed to vary through their thickness direction according to the power-law distribution of the volume fractions of the constituents. The von Kármán nonlinear strain–displacement relationship has been employed, and the shape functions of the approximations of the displacement fields constructed in terms of a set of mesh-free kernel particle functions. A SC nodal integration approach has been employed to obtain the bending stiffness, with the shear and membrane terms computed using a direct nodal integration technique to eliminate membrane and shear locking effects for very thin plates. The results obtained in these experiments have been compared

Fig. 12. The postbuckling response of aluminium/zirconia FGM plates (SSSS-01, SSSS02, a/h = 100) subject to edge compression.

Fig. 14. The postbuckling response of aluminium/zirconia FGM plates subject to edge compression and different temperature fields (SSSS-02, a/h = 20).

6. Conclusions

Y.Y. Lee et al. / Computer Methods in Applied Mechanics and Engineering 199 (2010) 1645–1653

with those obtained in the literature, and a good agreement is found. Numerical results have revealed that the FGM plates associated with a larger value of volume fraction exponent endure a larger postbuckling deformation than those with a smaller volume fraction exponent, and the existence of a temperature field reduces the buckling resistance of FGM plates and exacerbates the postbuckling deformation. Acknowledgement The work described in this paper was fully supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China [Project no. 9041239 (CityU 115607)]. References [1] G.N. Praveen, J.N. Reddy, Nonlinear transient thermoelastic analysis of functionally graded ceramic–metal plates, Int. J. Solids Struct. 33 (1998) 4457–4476. [2] N. Noda, Thermal stresses in functionally graded materials, J. Therm. Stress. 22 (1999) 477–512. [3] J. Yang, H.S. Shen, Vibration characteristics and transient responses of sheardeformable functionally graded plates in thermal environments, J. Sound Vibration 255 (2002) 579–602. [4] X.Q. He, T.Y. Ng, S. Sivashanker, K.M. Liew, Active control of FGM plates with integrated piezoelectric sensors and actuators, Int. J. Solids Struct. 38 (2001) 1641–1655. [5] K.M. Liew, X.Q. He, T.Y. Ng, S. Kitipornchai, Finite element piezothermoelasticity analysis and the active control of FGM plates with integrated piezoelectric sensors and actuators, Comput. Mech. 31 (2003) 350–358. [6] E. Efraim, M. Eisenberger, Exact vibration analysis of variable thickness thick annular isotropic and FGM plates, J. Sound Vibration 299 (2007) 720–738. [7] J.N. Reddy, Analysis of functionally graded plates, Int. J. Numer. Meth. Engrg. 47 (2000) 663–684. [8] S.S. Vel, R.C. Batra, Three-dimensional exact solution for the vibration of functionally graded rectangular plates, J. Sound Vibration 272 (2004) 703–730. [9] M.M. Najafizadeh, H.R. Heydari, Thermal buckling of functionally graded circular plates based on higher order shear deformation plate theory, Eur. J. Mech. A/Solids 23 (2004) 1085–1100. [10] K.M. Liew, J. Yang, S. Kitipornchai, Postbuckling of piezoelectric FGM plates subject to thermo-electro-mechanical loading, Int. J. Solids Struct. 40 (2003) 3869–3892. [11] H.S. Shen, Thermal postbuckling behavior of shear deformable FGM plates with temperature-dependent properties, Int. J. Mech. Sci. 49 (2007) 466–478. [12] J. Woo, S.A. Meguid, J.C. Stranart, K.M. Liew, Thermomechanical postbuckling analysis of moderately thick functionally graded plates and shallow shells, Int. J. Mech. Sci. 47 (2005) 1147–1171. [13] K.S. Na, J.H. Kim, Three-dimensional thermomechanical buckling analysis of functionally graded composite plates, Compos. Struct. 73 (2006) 413–422. [14] K.S. Na, J.H. Kim, Three-dimensional thermal buckling analysis of functionally graded materials, AIAA J. 43 (2005) 1605–1612. [15] K.S. Na, J.H. Kim, Thermal postbuckling investigation of functionally graded plates using 3-D finite element method, Finite Elem. Des. Anal. 42 (2006) 748–756. [16] C.Y. Chia, Nonlinear analysis of plates, McGraw-Hill, New York, 1980. [17] N. Yamaki, Postbuckling behavior of rectangular plates with small initial curvature loaded in edge compression, J. Appl. Mech. 26 (1959) 407–414. [18] H.S. Shen, Thermomechanical post-buckling analysis of imperfect laminated plates using a higher-order shear-deformation theory, Comput. Struct. 66 (1998) 395–409.

1653

[19] S. Wang, D.J. Dawe, Finite strip large deflection and post-overall-buckling analysis of diaphragm-supported plate structures, Comput. Struct. 61 (1996) 155–170. [20] B.Z. Huang, S.N. Atluri, A simple method to follow post-buckling paths in finite element analysis, Comput. Struct. 57 (1995) 477–489. [21] E. Madenci, A. Barut, Pre- and postbuckling responses of curved, thin, composite panels with cutouts under compression, Int. J. Numer. Meth. Engrg. 37 (1994) 1499–1510. [22] Y.H. Kim, A.K. Noor, Buckling and postbuckling of composite panels with cutouts subjected to combined loads, Finite Elem. Anal. Des. 22 (1996) 163–185. [23] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods for static and dynamic fracture, Int. J. Solids Struct. 32 (1995) 2547–2570. [24] P. Krysl, T. Belytschko, Analysis of thin plates by the element-free Galerkin method, Comput. Mech. 16 (1995) 1–10. [25] P. Krysl, T. Belytschko, Analysis of thin shells by the element-free Galerkin method, Int. J. Solids Struct. 33 (1996) 3057–3080. [26] S.N. Atluri, J.Y. Cho, H.G. Kim, Analysis of thin beams, using the meshless local Petrov–Galerkin method, with generalized moving least squares interpolations, Comput. Mech. 24 (1999) 334–347. [27] S. Li, W. Hao, W.K. Liu, Numerical simulations of large deformation of thin shell structures using meshfree method, Comput. Mech. 25 (2000) 102–116. [28] J.S. Chen, C. Pan, C.M.O.L. Roque, H.P. Wang, A Lagrangian reproducing kernel particle method for metal forming analysis, Comput. Mech. 22 (1998) 289–307. [29] K.M. Liew, L.X. Peng, S. Kitipornchai, Geometric non-linear analysis of folded plate structures by the spline strip kernel particle method, Int. J. Numer. Meth. Engrg. 71 (2007) 1102–1133. [30] X. Zhao, K.M. Liew, T.Y. Ng, Vibration analysis of laminated composite cylindrical panels via a mesh-free approach, Int. J. Solids Struct. 40 (2003) 161–180. [31] X. Zhao, K.M. Liew, T.Y. Ng, The element-free kp-Ritz method for free vibration analysis of conical shell panels, J. Sound Vibration 295 (2006) 906–922. [32] K.M. Liew, Y.G. Hu, X. Zhao, T.Y. Ng, Dynamic stability analysis of composite laminated cylindrical shells via the mesh-free kp-Ritz method, Comput. Meth. Appl. Mech. Engrg. 196 (2006) 147–160. [33] J.S. Chen, C.T. Wu, S. Yoon, Y. You, A stabilized conforming nodal integration for Galerkin mesh-free methods, Int. J. Numer. Meth. Engrg. 50 (2001) 435–466. [34] M.A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, John Wiley and Sons, Chichester, England, 2000. [35] J.N. Reddy, Mechanics of laminated composite plates and shells: theory and analysis, 2nd edCRC Press, Boca Raton, FL, 2004. [36] X. Zhao, Y.Y. Lee, K.M. Liew, Free vibration analysis of functionally graded plates using the element-free kp-Ritz method, J. Sound Vibration 319 (2009) 918–939. [37] J.S. Chen, C. Pan, C.T. Wu, W.K. Liu, Reproducing kernel particle methods for large deformation analysis of nonlinear structures, Comput. Meth. Appl. Mech. Engrg. 139 (1996) 237–262. [38] W.K. Liu, S. Jun, Y.F. Zhang, Reproducing kernel particle methods, Int. J. Numer. Meth. Engrg. 20 (1995) 1081–1106. [39] J. Dolbow, T. Belytschko, Nodal integration of element-free Galerkin method, Comput. Mech. 23 (1999) 219–230. [40] S. Beissel, snd Belytschko T, Nodal Integr. Element-free Galerkin Method Computer Methods in Applied Mechanics and Engineering 139 (1996) 49–74. [41] D.D. Wang, J.S. Chen, Locking-free stabilized conforming nodal integration for meshfree Mindlin–Reissner plate formulation, Comput. Meth. Appl. Mech. Engrg. 193 (2004) 1065–1083. [42] X. Zhao, G.R. Liu, K.Y. Dai, Z.H. Zhong, G.Y. Li, X. Han, A linearly conforming radial point interpolation method (LC-RPIM) for shells, Comput. Mech. 43 (2009) 321–441. [43] G.R. Liu, Y. Li, K.Y. Dai, M.T. Luan, W. Xue, A linearly conforming RPIM for solids mechanics problems, Int. J. Comput. Meth. 3 (2006) 401–428. [44] D.J. Dawe, S. Wang, Postbuckling analysis of thin rectangular laminated plates by spline FSM, Thin-Walled Struct. 30 (1998) 159–179.