Isogeometric analysis of functionally graded plates using higher-order shear deformation theory

Composites: Part B xxx (2013) xxx–xxx

Contents lists available at SciVerse ScienceDirect

Composites: Part B journal homepage: www.elsevier.com/locate/compositesb

Isogeometric analysis of functionally graded plates using higher-order shear deformation theory Loc V. Tran a, A.J.M. Ferreira b,c, H. Nguyen-Xuan a,d,⇑ a

Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam Universidade do Porto, Faculdade de Engenharia, Porto, Portugal c King Abdulaziz University, Jeddah, Saudi Arabia d Department of Mechanics, Faculty of Mathematics & Computer Science, University of Science, , Ho Chi Minh City, Viet Nam b

a r t i c l e

i n f o

Article history: Received 7 February 2013 Accepted 23 February 2013 Available online xxxx Keywords: A. Laminates B. Buckling C. Computational modelling

a b s t r a c t This paper presents a novel and effective formulation based on isogeometric approach (IGA) and higherorder deformation plate theory (HSDT) to study the behavior of functionally graded material (FGM) plates. HSDT model using C1 continuous element is able to improve the accuracy of solution and describe exactly the shear stress distribution without shear correction factors. IGA utilizes the non-uniform rational B-spline (NURBS) functions which allow us to achieve easily the smoothness with arbitrary continuity order. The present method hence fulfills the C1 – requirement of HSDT model. The effective material properties of the FGM plates, which property varies only through the thickness of plate, are calculated using the rule of mixture and the Mori–Tanaka homogenization technique. The static, dynamic and buckling analysis of rectangular and circular plates is investigated for different boundary conditions. Numerical results show high effectiveness of the present formulation. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Laminated composite material is made by stacking several lamina layers together to achieve desired properties e.g. high stiffness and strength-to-weight ratios, long fatigue life, wear resistance, etc. For example, a ceramic layer with the low thermal conductivity when bonded to a metal structure with high ductility can resist high thermal environment or extremely large temperature gradients. Therefore, they are more suitable to use in aerospace structure applications and nuclear plants. However, the sudden changed material at the interface causes the large inter–laminar stresses leading to delamination [1]. To overcome this shortcoming, in 1984, a group of scientists in Sendai – Japan proposed a new smart material so–called a functionally graded material (FGM) [2–5] with continuous change of material properties along certain dimensions of the structure. This paper follows the FGM plate – a mixture of two distinct material phases: ceramic and metal with the variation of the volume fraction according to power index through the plate thickness. The effective material properties such as Young’s modulus, Poisson’s ratio, density, thermal conductivity, and thermal expansion can be determined by several rules such as the mixture [6], the three-phase model by Frohlich and

⇑ Corresponding author at: Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam. E-mail address: [email protected] (H. Nguyen-Xuan).

Sack [7], the self-consistent scheme [8], the Mori–Tanaka technique [9], and the mean field approach [10]. The theoretical frame of FGM plates is commonly relied on popular plate theories such as the Classical Plate Theory (CPT), First Order Shear Deformation Theory (FSDT) and High Order Shear Deformation Theory (HSDT). The CPT relied on the Kirchhoff–Love assumptions ignores the transverse shear deformation. Hence, it cannot produce accuracy results for moderate and thick plates (with thickness/span > 1/10). The FSDT [11–13], which takes into account the effects of shear deformation, was therefore developed. The FSDT is simple to implement and applied for both thick and thin FGM plates. However, the accuracy of solutions will be strongly dependent on the shear correction factors in which their values are quite dispersed through each problem. In addition, with constant shear strain assumption, shear stress obtained from FSDT distributes wrongly and is not satisfactory the traction free boundary conditions on the plate’s surfaces. The HSDT [6,14,15,52] that includes higher-order terms in the approximation of the displacement field has then been developed. This theory disregards shear correction factors and yields more accurate and stable solutions (e.g., inter-laminar stresses and displacements [16,17]) than the FSDT ones. It is worth mentioning that the HSDT requires C1-continuity of generalized displacements leading to the second-order derivative of the stiffness formulation and hence they cause some obstacles in standard finite formulations based on the C0-continuity. In this paper, we show that such a C1-HSDT formulation will be easily achieved using a NURBS-based isogeometric approach.

1359-8368/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compositesb.2013.02.045

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

2

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

Isogeometric approach, a new computational method, has been recently proposed by Hughes et al. [18] to closely link the gap between Computer Aided Design (CAD) and Finite Element Analysis (FEA). The crucial idea is that the IGA uses the same non-uniform rational B-spline (NURBS) functions in describing the geometry of problem in CAD and constructing finite approximation for analysis. The exact geometry is therefore maintained at the coarsest level of discretization and the re-meshing is performed on this coarsest level without any communication with CAD geometry. Furthermore, NURBS functions provide a flexible way to make refinement, derefinement, and degree elevation [19]. They enable us to achieve easily the smoothness of arbitrary continuity order in comparison with the traditional FEM. Hence, IGA naturally verifies the C1-continuity of plates based on the HSDT assumption. Isogeometric analysis has been widely applied to various practical problems such as linear and non-linear elasticity and plasticity behavior [20], structural vibrations [21], the Reissner–Mindlin plate and shell [22,23], Kirchhoff–Love shell [24,25], structural shape optimization [26], and further improved NURBS approaches [27,28], just mention to a few. This paper presents a NURBS-based isogeometric approach for static, dynamic and buckling analysis of the FGM plates using HSDT model. The material property changing continuously through plate thickness is homogenized by the rule of mixture and the Mori–Tanaka homogenization technique. IGA utilizes the non-uniform rational B-spline (NURBS) functions which achieve easily the smoothness with arbitrary continuity order. The method satisfies naturally the C1 continuous requirement of this plate theory. Several numerical examples are illustrated to demonstrate high accuracy of the present method in comparison with other published methods. The paper is outlined as follows. Next section describes a FGM plate model. In Section 3, we introduce a novel plate formulation based on IGA and HSDT. The numerical results and discussions are provided in Section 4. Finally, the article is closed with some concluding remarks.

where subscripts m and c refer to the metal and ceramic constituents, respectively. Eq. (1) implies that the volume fraction varies through the thickness based on the power index n. The effective material properties according to the rule of mixture are the given by

PðzÞ ¼ ðPc  Pm ÞV c ðzÞ þ Pm

ð2Þ

where Pc, Pm denote the material properties of the ceramic and the metal, respectively, including the Young’s modulus E, Poisson’s ratio m, density q. However, the rule of mixture does not consider the interactions among the constituents [29]. So, the Mori–Tanaka technique [9] is then developed to take into account these interactions with the effective bulk and shear modulus defined from following relations:

Ke  Km Vc ¼ c K m K c  K m 1 þ V m K Kþ4=3 l m

m

ð3Þ

le  lm Vc ¼ lc  lm 1 þ V m llcmþflm1

ð9K m þ8lm Þ where f1 ¼ lm6ðK . And the effective values of Young’s modulus m þ2lm Þ E and Poisson’s ratio m are given by

E¼

9K e le ; 3K e þ le

m¼

3K e  2le 2ð3K e þ le Þ

ð4Þ

Fig. 2 illustrates comparison of the effective Young’s modulus of Al/ZrO2 FGM plate calculated by the rule of mixture and the Mori– Tanaka scheme via the power index n. Note that with homogeneous material, two models produce the same values. As material is inhomogeneous, the effective property through the thickness of the former is higher than the latter one. 2.2. Governing equation for third-order shear deformation theory model Based on third-order shear deformation theory [30], the displacements of an arbitrary point in the plate are defined as

2. Governing equations and weak forms for FGM plates

uðx; y; zÞ ¼ u0 þ zbx þ cz3 ðbx þ w;x Þ

2.1. Functionally graded material

v ðx; y; zÞ ¼ v 0 þ zby þ cz3 ðby þ w;y Þ;

Functionally graded material is a composite material which is created by mixing two distinct material phases. Two mixed materials are often ceramic at the top and metal at the bottom as shown in Fig. 1. In our work, two homogenous models have been used to estimate the effective properties of the FGM include the rule of mixture [6] and the Mori–Tanaka technique [9]. Herein, we adopt the volume fraction of the ceramic and metal phase by the power–law exponent function

wðx; yÞ ¼ w0

 n 1 z V c ðzÞ ¼ ; þ 2 h

Vm ¼ 1  Vc

ð1Þ

  h h 6z6 2 2

ð5Þ

where c = 4/3h2 and the variables u0 = {u0 v0}T, w0 and b = {bx by}T are the membrane displacements, the deflection of the mid-plane and the rotations in the x–z, y–z planes, respectively. The relationship between strains and displacements is described by

e ¼ ½exx eyy cxy T ¼ e0 þ zj1 þ z3 j2 c ¼ ½cxz cyz T ¼ es þ z2 js

ð6Þ

where

2

e0 ¼ 6 4

u0;x

3

2

bx;x

3

2

bx;x þ w0;xx

3

v 0;y 75; j1 ¼ 64 by;y 75; j2 ¼ c64 by;y þ w0;yy 75 bx;y þ by;x bx;y þ by;x þ 2w0;xy u0;y þ v 0;x ð7Þ

and

"

es ¼

bx þ w0;x by þ w0;y

#

" ; js ¼ 3c

bx þ w0;x by þ w0;y

# ð8Þ

A weak form of the static model for the plates under transverse loading q0 can be briefly expressed as [12]: Fig. 1. The functionally graded plate model.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

3

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

0.5 n=10

0.4 0.3

n=5

0.2 n=3

0.1

n=1

0 n=0.5

−0.1

n=0.3

−0.2 −0.3

n=0.1

−0.4

n=0

−0.5 1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

Fig. 2. The effective modulus of Al/ZrO2 FGM plate computed by the rule of mixture (in solid line) and the Mori–Tanaka (in dash dot line).

Z

deT DedX þ

Z

X

dcT Ds cdX ¼

X

Z

dw0 q0 dX

ð9Þ

X

where

2

A

B

6 D ¼ 4B E

3

E

"

7 D F 5; F H

DS ¼

AS

BS

BS

DS

ð16Þ

# ð10Þ

A¼ E¼

Z

h=2

Q dz;

B¼

Z

h=2 h=2

F¼

h=2

As ¼

Z

h=2

Q zdz;

D¼

Z

h=2

Q z3 dz; h=2

Bs ¼

Gdz;

For the buckling analysis, a weak form of the plate under the inplane forces can be expressed as:

Z

in which

Z

8 9 8 9 8 8 9 9 > > > > < w0 > < u1 > < u0 > < v0 > = = = = ~ ¼ u2 ; u1 ¼ by bx u ; u2 ¼ ; u3 ¼ 0 > > > > > > : > : > : : ; ; ; ; by þ w0;y bx þ w0;x 0 u3

h=2

Z

Z

Q z2 dz;

where "

h=2

h=2

Q z4 dz;

H¼

h=2 h=2

Ds ¼

Gz2 dz;

Z

h=2

Z

X

h=2

Z

deT DedX þ

h=2

Q z6 dz

ð11Þ

h=2

N0 ¼

dcT Ds cdX þ

Z

X

rT = [@/@x N 0x N 0xy

N0xy N0y

#

rT dw0 N0 rw0 dX ¼ 0

ð17Þ

X

o/oy]T

is

the

gradient

operator

and

is a matrix related to the pre-buckling loads.

h=2

Gz4 dz

3. A novel plate formulation based on isogeometric approach and the C1-HDST type

h=2

where the material matrices are given as

2 3 1 mðzÞ 0 EðzÞ 6 7 Q¼ 1 0 4 mðzÞ 5 1  m2 ðzÞ 0 0 ð1  mðzÞÞ=2   1 0 EðzÞ G¼ 2ð1 þ mðzÞÞ 0 1

3.1. Knot-vector and basis function

ð12Þ

and the Poisson’s ratio m(z) and Young modulus E(z) are the effective material properties and vary along the thickness direction. For the free vibration analysis of the plates, weak form can be derived from the following dynamic equation

Z

deT DedX þ

X

Z

dcT Ds cdX ¼

X

Z

€~ dX ~ T mu du

ð13Þ

X

where m the mass matrix is calculated according to consistent form

2

I0 6 m¼40

0 I0

3 2 I1 0 7 6 0 5 where I0 ¼ 4 I2

I2 I3

0

0

I0

cI5

ðI1 ; I2 ; I3 ; I4 ; I5 ; I7 Þ ¼

cI4 Z

qðzÞð1; z; z2 ; z3 ; z4 ; z7 Þdz

 Ni;0 ðnÞ ¼

Ni;p ðnÞ ¼

ð15Þ

1 if ni 6 n < niþ1 0 otherwise

 ð18Þ

The basis functions are defined by the following recursion formula p P 1

ð14Þ

h=2

h=2

and

3 cI4 7 cI5 5 2 c I7

A knot vector N = {n1, n2, . . . , nn+p+1} is a non-decreasing sequence of parameter values ni ; i ¼ 1; . . . ; n þ p, where ni 2 R called ith knot lies in the parametric space, p is the order of the B-spline and n is number of basis functions needed to describe the B-spline curve. If the first and the last knots are repeated p + 1 times, the knot vector is called open knot. A B-spline basis function is C1 continuous inside a knot span and Cp1 continuous at a single knot. The B-spline basis functions Ni,p(n) are defined recursively on the corresponding knot vector start with order p = 0 as follows

niþpþ1  n n  ni Ni;p1 ðnÞ þ Niþ1;p1 ðnÞ niþp  ni niþpþ1  niþ1

ð19Þ

It is well known that for p = 0 or 1 the basis functions of isogeometric analysis are identical to those of standard piecewise constant and linear finite elements, respectively. However, they are different as p P 2. In this study, we consider the basis functions

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

4

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

1

1

0.9

0.9

0.8

0.8

0.7

0.7

p=1

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

p=2

0.5

1

0

0

1

1

0.9

0.9

0.8

1

0.8

p=3

0.7

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1 0

1/2

p=4

0.7

0.6

0

1/2

1

0

0

1/2

1

Fig. 3. Some B-spline basis functions: linear, quadratic, cubic and quartic functions.

with p P 2. Fig. 3 illustrates a set of one-dimensional linear, quadratic, cubic and quartic B-spline basis functions according to open uniform knot vectors N = {0, 0, 0.5, 1, 1}, N = {0, 0, 0, 0.5, 1, 1, 1}, N = {0, 0, 0, 0, 0.5, 1, 1, 1, 1}, N = {0, 0, 0, 0, 0, 0.5, 1, 1, 1, 1, 1}, respectively. It is clear that the basis functions are Cp1 continuous. Thus, as p P 2, the present approach always satisfies C1-requirement in approximate formulations based on the third shear deformation plate theory proposed by Reddy [30].

where NbI ðn; gÞ ¼ Ni;p ðnÞMj;q ðgÞ is the shape function associated with node I. To present exactly curved geometries (e.g., circles, cylinders, spheres, etc.) non-uniform rational B-splines (NURBS) are used. Be different from B-spline, each control point of NURBS has additional value called an individual weight wI. Then the NURBS surface can be expressed as

Sðn; gÞ ¼

mn X

NI ðn; gÞPI

I¼1

3.2. B-spline curves and surface

n X Ni;p ðnÞPi

ð23Þ

I¼1

The B-spline curve is defined as

CðnÞ ¼

N b wI with NI ðn; gÞ ¼ mn I X b NI ðn; gÞwI

ð20Þ

It notes that B-spline function is only the special case of the NURBS function when the individual weight of control points is constant. An exact cubic NURBS surface with its control points in red1 color is illustrated in Fig. 4.

i¼1

where Pi are the control points and Ni,p(n), which is the pth degree B-spline basis function, is defined on the open knot vector. The B-spline surfaces are defined by the tensor product of basis functions in two parametric dimensions n and g with two knot vectors N = {n1, n2, . . . , nn+p+1} and H = {g1, g2, . . . , gm+q+1} are expressed as follows

3.3. NURBS-based novel plate formulation Using the NURBS basis functions above, the displacement field u of the plate is approximated as

uh ðn; gÞ ¼

mn X

NI ðn; gÞqI

ð24Þ

I¼1

Sðn; gÞ ¼

n X m X

Ni;p ðnÞMj;q ðgÞPi;j

ð21Þ

i¼1 j¼1

where Ni,p(n) and Mj,q(g) are the B-spline basis functions defined on the knot vectors over an m  n net of control points Pi,j. Similarly notations used in finite elements, Eq. (21) can be rewritten as

Sðn; gÞ ¼

m n X

N bI ðn; gÞPI

I¼1

ð22Þ

where qI = [u0I v0I w0I bxI byI]T is the vector of nodal degrees of freedom associated with the control point I. Substituting Eq. (24) into Eqs. (7) and (8), the in-plane and shear strains can be rewritten as:

1 For interpretation of color in Fig. 4, the reader is referred to the web version of this article.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

5

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

Substituting Eq. (25) into Eqs. (9), (13), and (17), the formulations of static, free vibration and buckling problem are rewritten in the following form

Kd ¼ Fl

ð27Þ

ðK  x2 MÞd ¼ 0

ð28Þ

ðK  kcr Kg Þd ¼ 0

ð29Þ

where the global stiffness matrix K is given by

8 m 9T 2 38 m 9 ( )T " Z > = A B E > = As Bs0 6 7 K¼ Bb1 4 B D F 5 Bb1 þ s1 > Bs B X> : b2 > : b2 > ; ; E F B B B

Bs

#(

Ds

Bs0 Bs1

) dX

ð30Þ Fig. 4. Physical mesh and control mesh with cubic non-uniform rational B-spline surface.

and the load vector is computed by

Fl ¼

Z

q0 NdX

ð31Þ

X

the global mass matrix M is expressed as

Table 1 Material properties. Al

SiC

ZrO2  1

ZrO2  2

Al2O3

M¼

Z

e T m Nd e X N

ð32Þ

X

E (GPa)

70 0.3 2707

m q (kg/m3)



T

eT0 jT1 jT2 eTs jTs ¼

427 0.17 – mn Xh

200 0.3 5700

151 0.3 3000

380 0.3 3800

T b1 T b2 T s0 T s1 T ðBm I Þ ðBI Þ ðBI Þ ðBI Þ ðBI Þ

iT

qI

I¼1

ð25Þ in which

2

NI;x

0

0 0 0

3

2

0 0 0

6 7 6 NI;y 0 0 0 5; Bb1 Bm I ¼ 4 0 I ¼ 40 NI;y NI;x 0 0 0 0 2 3 0 0 NI;xx NI;x 0 6 7 0 NI;yy 0 NI;y 5; Bb2 I ¼ c4 0 0 0 2NI;xy NI;y NI;x    0 0 NI;x NI 0 0 s0 ; Bs1 BI ¼ I ¼ 3c 0 0 NI;y 0 NI 0

0 0

0

0 0 NI;y

8 9 2 NI 0 0 0 > < N1 > = e ¼ N2 ; N1 ¼ 6 N 4 0 NI 0 0 > : > ; 0 0 NI 0 N3 2 3 2 0 0 0 NI 0 0 6 7 6 N2 ¼ 4 0 0 0 0 NI 5; N3 ¼ 4 0 0 0 0

3

0

NI;x

where

0

0

3

7 0 5; 0 NI

0 NI;y

0

7 NI 5

0

0

the geometric stiffness matrix is

NI;x

Kg ¼

0

ðBg ÞT N0 Bg dX

0

3

NI;x

0

0 0

7 NI;y 5;

Z

0

ð33Þ

ð34Þ

X

where

0

NI;x

NI

0

0 NI;y

0

NI



BgI ¼ ð26Þ



0 0

NI;x

0 0

0 0 NI;y

0 0

 ð35Þ

in which x, kcr 2 R are the natural frequency and the critical buckling value, respectively. 1 0.9 0.8 0.7

y

0.6 0.5 0.4 0.3 0.2 0.1 0

(a)

0

0.2

0.4

x

0.6

0.8

1

(b)

Fig. 5. Square plate: (a) model; (b) meshing of 11  11 cubic B-spline elements.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

6

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

1.53

1.528

Transerve displacement

1.5275

Deflection

1.467

1.52

Ref IGA p = 2 IGA p = 3

1.527

HSDT

1.5265

1.51

1.4665 Ref

1.5

IGA p = 2 IGA p = 3

1.466

1.49

FSDT 1.48

1.4655 5

10

15

20

25

1.47

5

10

15

20

1.46 5

25

10

15

20

25

Number element per edge

Number element per edge

(a) Power index n = 1 2.35

Transerve displacement

2.3215

Deflection

2.321 Ref IGA p = 2 IGA p = 3

2.3205

HSDT

2.32

2.3

2.16

2.25

Ref IGA p = 2 IGA p = 3

2.1595

FSDT

2.2 2.159 5

10

15

20

25

2.3195 5

10

15

20

25

2.15 5

10

Number element per edge

15

20

25

Number element per edge

(a) Power index n = 6 Fig. 6. Comparison of present result with analytical solution of Vel and Batra according to power index n = 1 and 6, respectively.

0.5

0.5

0.4

0.4

n = 0, HSDT FSDT n = 2, HSDT FSDT

0.3 0.2

0.3

0.1

0.1

0

0

−0.1

−0.1

−0.2

−0.2

−0.3

−0.3

−0.4

−0.4

−0.5 −0.25

−0.2

−0.15

−0.1

−0.05

0

0.05

0.1

n = 0, HSDT FSDT n = 2, HSDT FSDT

0.2

0.15

−0.5 −0.4

−0.35

−0.3

−0.25

−0.2

−0.15

−0.1

−0.05

0

Fig. 7. The stresses through thickness of Al/ZrO2  1 FG plate under uniform loads a/h = 5 and clamped edges based on HSDT and FSDT.

It is observed from Eq. (30) that the shear correction factors are no longer required and the stiffness formulation involves additional terms based on the third shear deformation plate theory

[30]. Herein Bb2 I contains the second-order derivative of the shape function. Hence, it requires C1-continuity in approximate formulations. It is worth mentioning that the construction of a C1 continu-

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

7

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

1.5

CCCC: IGA MLPG SCSC: IGA MLPG SSSS: IGA MLPG SFSF: IGA MLPG

Deflection

1

0.5

0

0 0.5 1

2

metal

Power index Fig. 8. The normalized deflection of Al/ZrO2  1 FGM plate via power indexes and boundary conditions.

Fig. 9. Deflection of Al/ZrO2  1 FG plate via boundary conditions: (a) SFSF; (b) SSSS; (c) CSCS; (d) CCCC.

ous finite element formulation having five-degrees-of-freedom (DOF) per node is more complicated and difficult. A C0 continuous element [16,32] with seven-degrees-of-freedom (DOF) was then proposed. However, the increasing of the computational cost and the accuracy of solution are still questionable. It is now interesting that our present formulation based on isogeometric analysis satisfies naturally from the theoretical/mechanical viewpoint of FGM plates [33]. In our work, the basis functions are Cp1 continuous.

As a result, as p P 2, the present approach always satisfies C1requirement in approximate formulations based on the proposed third shear deformation plate theory. 4. Results and discussion In this section, using IGA, we study the behavior of the FGM plates with rectangular and circular shapes. Herein, ceramic–me-

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

8

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

0.5

0.5 IGA (p=3) MLPG

0.4 0.3

0.3

0.2

0.2

0.1

0.1

0

0

Ceramic or metal

−0.1

−0.2

−0.3

−0.3

−0.4

−0.4 −0.15

−0.1

−0.05

0

0.05

0.1

0.15

0.5

−0.5 −0.25 −0.2 −0.15 −0.1 −0.05

0

0.05

0.1

0.15

0.5 IGA (p=3) MLPG

0.4

0.3

0.2

0.2

0.1

0.1

0

IGA (p=3) MLPG

0.4

0.3

−0.1

n = 0.5

−0.1

−0.2

−0.5 −0.2

IGA (p=3) MLPG

0.4

0

n=1

−0.1

−0.2

−0.2

−0.3

−0.3

−0.4

−0.4

−0.5 −0.25 −0.2 −0.15 −0.1 −0.05

0

0.05

0.1

0.15

n=2

−0.5 −0.25 −0.2 −0.15 −0.1 −0.05

0

0.05

0.1

0.15

Fig. 10. Comparison of the axial stress through thickness of Al/ZrO2  1 FG plate under uniform load, a/h = 5 and clamped edges.

tal functionally graded plates, which material properties are given in Table 1, are considered. Numerical results are obtained using full (p + 1)  (q + 1) Gauss points. 4.1. Static analysis 4.1.1. Comparison study Let us consider the simply supported Al/SiC square FGM plate shown in Fig. 5a which properties are given in Table 1. The plate



 is subjected to a sinusoidal pressure defined as q0 sin pax sin pay at the top surface. A convergent study of transverse displacement by quadratic (p = 2) and cubic (p = 3) elements is depicted in Fig. 6 according to n = 1 and 6, respectively. It is observed that the IGA gains the supper convergence. The discrepancy between meshing of 5  5 and 25  25 is lower than 0.1%. The results using HSDT model are better than that of FSDT one and match very well with the reference solution using 3D deformation model by Vel and Batra [1]. Here, just using 11  11 cubic NURBS elements as shown in Fig. 5b, IGA using HSDT produces an ultra-accurate solution that is very close to the exact one with very small error around 0.5‰. In addition, the cubic element with the higher order approximation obtains the better results compared to quadratic one. Next, we illustrate the performance of the present formulation using a specified degree of p = 3 (or cubic NURBS element). Both FSDT and HSDT theories are used to investigate the stress distributions through the thickness of the Al/ZrO2  1 FGM plate with a/ h = 5. Fig. 7a plots the axial stress at the central point of the clamped plate according to n = 0, 2, respectively. It is revealed that

Table 2  c ¼ 100wc =h of square Al/ZrO2  2 plate subNon-dimensional center deflection w jected to a uniform load q0. a/h

Method

20

Meshless method MLPG IGA (p = 2) IGA (p = 3) IGA (p = 4)

5

Meshless method MLPG IGA (p = 2) IGA (p = 3) IGA (p = 4)

N 0

0.5

1

2

1

2.08 (2.08) 2.118 2.0827 (2.0827) 2.0831 (2.0831) 2.0831 (2.0831)

2.79 2.65 – 2.7919 (2.6611) 2.7924 (2.6616) 2.7924 (2.6616)

3.09 2.97 3.150 3.0916 (2.9778) 3.0922 (2.9783) 3.0922 (2.9783)

3.33 3.24 3.395 3.3392 (3.2512) 3.3398 (3.2518) 3.3398 (3.2518)

4.48 (4.48) 4.580 4.4928 (4.4928) 4.4935 (4.4935) 4.4935 (4.4935)

2.477 (2.477) 2.436 2.4818 (2.4818) 2.4819 (2.4819) 2.4819 (2.4819)

3.2930 (3.135) – 3.3008 (3.1424) 3.3009 (3.1425) 3.3009 (3.1425)

3.6660 (3.515) 3.634 3.6739 (3.5231) 3.6741 (3.5233) 3.6741 (3.5233)

4.0090 (3.883) 3.976 4.0160 (3.8903) 4.0162 (3.8905) 4.0162 (3.8905)

5.343 (5.343) 5.253 5.3536 (5.3536) 5.3538 (5.3538) 5.3538 (5.3538)

Results according to the rule of mixture are in parentheses.

with homogeneous material, axial stress is symmetric through the middle plan. In inhomogeneous case, the stress is distributed in tendency in which the magnitude of compressive stress at the top is greater than the tensile one at the bottom. And we can see there are the overlapping results between FSDT and HSDT model.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

9

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

0.8

(a)

(b)

4

0.6 0.4

7

1

0.2 0

8

5

2

−0.2

3

9

−0.4 −0.6

6 −0.8 −0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

Fig. 11. Circular plate: (a) geometry; (b) the coarsest mesh with only one quadratic element corresponding to nine control points in red color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3 The coordinates and weight values of control points of a circular plate R = 05. I xi yi wi

1

2

3

4

5

6

pffiffiffi  2=4 pffiffiffi 2=4 1

pffiffiffi  2=2 0 pffiffiffi 2=2

pffiffiffi  2=4 pffiffiffi  2=4 1

0 pffiffiffi 2=2 pffiffiffi 2=2

0

0

0 1

pffiffiffi  2=2 pffiffiffi 2=2

7 pffiffiffi 2=4 pffiffiffi 2=4 1

8 pffiffiffi 2=2 0 pffiffiffi 2=2

9 pffiffiffi 2=4 pffiffiffi  2=4 1

Nevertheless, this observation is no longer true in case of shear xz ða=2; 0Þ plotted in Fig. 7b. Applying FSDT, the shear stress stress s is wrong distribution with constant (n = 0) and negative values at two free surfaces. While HSDT using third shear deformation theory satisfies the condition of parabolically curved shear stress distribution with traction free boundary at top and bottom surfaces of the plates [31]. 4.1.2. Behavior of square plate under various boundary conditions The effective of boundary conditions on the normalized central  c ¼ 100wc Em h3 =f12ð1  m2 Þp0 a4 g of Al/ZrO2  1 plate is deflection w

investigated in Fig. 8. It is seen that the cubic element combined with HSDT approach gives the results that are very close to MLPG method by Gilhooley et al. [34] for all boundary conditions. Moreover, when the boundary condition changes from CCCC to SCSC, SSSS and SFSF, the structural stiffness is reduced, the magnitude of deflection is thus increased, respectively. The shapes of transverse displacement according the various boundary conditions are illustrated in Fig. 9. Fig. 10 reveals the normalized axial stresses r x through the plate thickness with volume exponent n = 0, 0.5, 1, 2, respectively. It can be seen that the obtained results match well with those of MLPG approach presented by Gilhooley et al. [34].

4.1.3. Effect of the types of homogeneous scheme Table 2 summarizes IGA results in comparison with those of meshless method [14] and MLPG [29] using the same HSDT for thick and moderate Al/ZrO2  2 plates. In this example, both homogeneous schemes: the Mori–Tanaka homogenization method and the rule of mixture are used. The obtained non-dimensional deflections show some following remarks: (1) with meshing of

0.6

0.5 0.4

0.4

0.3 0.2

0.2

0

y

y

0.1 0

−0.1 −0.2

−0.2 −0.3

−0.4

−0.4 −0.6 −0.6

−0.4

−0.2

0

0.2

0.4

0.6

−0.5 −0.5

0

x

x

(a)

(b)

0.5

Fig. 12. The meshing of circle: (a) a cubic element; (b) 13  13 cubic elements.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

10

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

Fig. 13. Deflection of circular FGM plate with n = 2, h/R = 0.1: (a) clamped boundary; (b) roller supported boundary.

Table 4  c ¼ 100wc =h of circular Al/ZrO2  2 plates subjected to a uniform load q0. Non-dimensional center deflection w Boundary condition

Roller supported

Clamped

h/R

Method or author

Power index n 0

2

4

8

10

50

100

104

0.05

Li Reddy Present

10.391 10.396 10.2203

5.709 5.714 5.6100

5.219 5.223 5.1299

4.809 4.812 4.7299

4.70 4.704 4.6239

4.256 4.258 4.1881

4.187 4.189 4.1198

4.118 4.118 4.0495

0.1

Li Reddy Present

10.46 10.481 10.3440

5.738 5.756 5.6742

5.245 5.261 5.1879

4.835 4.848 4.7838

4.727 4.739 4.6769

4.283 4.293 4.2377

4.215 4.223 4.1690

4.146 4.151 4.0970

0.15

Li Reddy Present

10.575 10.623 10.4962

5.786 5.826 5.7494

5.289 5.325 5.2556

4.879 4.909 4.8478

4.771 4.799 4.7401

4.329 4.349 4.2987

4.26 4.28 4.2297

4.191 4.208 4.1572

0.2

Li Reddy Present

10.736 10.822 10.6973

5.853 5.925 5.8475

5.351 5.414 5.3439

4.941 4.993 4.9315

4.833 4.882 4.8230

4.392 4.429 4.3793

4.324 4.359 4.3098

4.255 4.286 4.2369

0.05

Li Reddy Present

2.561 2.554 2.5480

1.405 1.402 1.3990

1.284 1.282 1.2786

1.184 1.181 1.1785

1.157 1.155 1.1520

1.049 1.046 1.0435

1.032 1.029 1.0267

1.015 1.011 1.0092

0.1

Li Reddy Present

2.667 2.639 2.6297

1.456 1.444 1.4386

1.329 1.32 1.3143

1.227 1.217 1.2123

1.201 1.19 1.1855

1.091 1.08 1.0762

1.074 1.063 1.0592

1.057 1.045 1.0415

0.15

Li Reddy Present

2.844 2.781 2.7652

1.54 1.515 1.5043

1.405 1.384 1.3733

1.30 1.278 1.2685

1.273 1.25 1.2412

1.162 1.137 1.1304

1.145 1.119 1.1132

1.127 1.101 1.0952

0.2

Li Reddy Present

3.093 2.979 2.9541

1.658 1.613 1.5958

1.511 1.473 1.4557

1.402 1.362 1.3467

1.375 1.333 1.3187

1.262 1.216 1.2060

1.244 1.199 1.1884

1.226 1.18 1.1700

15  15 elements, the obtained results using cubic and quartic elements are almost the same, which offers to use only cubic element in what follows; (2) the present results are excellent agreement with those by meshless method [14] using 19  19 number of collocated points for both homogeneous schemes; (3) Fig. 2 shows that the rule of mixtures holds the larger effective Young modulus through the plate thickness than that of the Mori–Tanaka scheme. So, the transverse displacement by former model is slightly smaller compared to these of the latter one because of stiffer stiffness matrix. 4.1.4. Behavior of circular plate under uniform load In this section, let study behavior of FGM circular plate with h thickness and radius R = 0.5 (see in Fig. 11a) which are made of titanium/zirconium [35]. From the tensor product N  H with two knot vectors N = {0, 0, 0, 1, 1, 1} and H = {0, 0, 0, 1, 1, 1}, the geometry of circular plate is described exactly in Fig. 11b with only one quadratic element corresponding to nine control points which

coordinate is given in Table 3. By repeating the first and the last values in knot vectors one time, the polynomial order of basis functions is increased from p = 2 to p = 3. As a result, the exact geometry preserved with only one cubic element is plotted in Fig. 12a. Then applying the knot insertion technique, we can obtain the meshing of 13  13 cubic NURBS elements as depicted in Fig. 12b. Based on the rule of mixture, the effective modulus is estimated following Ref. [36] by Reddy et al.

E ¼ ðEm  Ec Þð0:5  z=hÞn þ Ec

ð36Þ

Eq. (36) shows that as n = 0, the plate is fully metal and as n ? 1, the homogeneous ceramic plate is retrieved. Under the uniform pressure at the top of plate, the deflection shapes of the plates are plotted in Fig. 13 with the normalized maximum values at the center summarized in Table 4. It is observed that the obtained results are in good agreement with that of literatures by Reddy et al. [36] and elasticity solution by Li et al. [35]. As boundary condition is simply soft supported (just w = 0 at boundary) the difference be-

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

11

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

0.5

0.5 Metal n=2 n=4 n = 10 Ceramic

0.4 0.3

0.4 0.3

0.2

0.2

0.1

0.1

0

0

−0.1

−0.1

−0.2

−0.2

−0.3

−0.3

(a)

−0.4 −0.5

−0.1

−0.4

−0.05

0

0.05

0.1

Metal n=2 n=4 n = 10 Ceramic

(b)

−0.5 −0.2

0.15

−0.15

−0.1

−0.05

0

Fig. 14. The stresses of circular FGM plate with clamped edge and R/h = 10 (a): axial stress; (b) shear stress.

0.5

0.5 Metal n=2 n=4 n = 10 Ceramic

0.4 0.3

0.4 0.3

0.2

0.2

0.1

0.1

0

0

−0.1

−0.1

−0.2

−0.2

−0.3

−0.3

−0.4

(a)

−0.5 −0.4

−0.3

−0.4 −0.2

−0.1

0

0.1

0.2

0.3

−0.5 −0.2

Metal n=2 n=4 n = 10 Ceramic

(b) −0.15

−0.1

−0.05

0

Fig. 15. The stresses of circular FGM plate SSSS R/h = 10 (a): axial stress at center; (b) shear stress at r = R/2.

Table 5 The comparison of natural frequency of simply supported Al/ZrO2 FG plate with a/h = 5. n

IGA 55

77

11  11

15  15

FSDT

1 2 5

0.2185 0.2190 0.2216

0.2185 0.2190 0.2216

0.2185 0.2190 0.2216

0.2185 0.2190 0.2216

0.2216 0.2230 0.2264

tween the present method and the elasticity solution is very small and reduces as a plate becomes thicker. In fact, it decreases from 1.6% to 0.5% as thickness to radius ratio increases in range from 0.05 to 0.2. In case of clamped plate, the IGA result is closer to that of Reddy et al. [36] than the elasticity solution with very small difference around 0.8%. Figs. 14 and 15 reveal the behavior of clamped FG plate under  r ¼ rr ð0; 0Þ=q0 and the uniform load consist of the radial stress r rz ¼ srz ðR=2; 0Þ=q0 under clamped and roller supshear stress s ported condition, respectively. It is also seen that, with isotropic plates, the linear axial stress holds zero at the mid-plane because of the symmetry. However, in the FG plates, because the components of the material particle are arranged in tendency of increasing the stiffness according to increment of z value, the axial stress is redistributed with increase of compression at the top (ceramic surface) and decrease of tension at the bottom (metal surface),

Meshless method [15]

MLPG [29]

Exact [1]

0.2188 0.2188 0.2215

0.2152 0.2153 0.2194

0.2192 0.2197 0.2225

respectively. In addition, the shear stress is parabolically distributed with the peak point where the maximum value is archived also tends to move to the ceramic surface.

4.2. Vibration analysis 4.2.1. Square plate First, let consider a simply supported Al/ZrO2  1 plate which is homogenized follow to Mori–Tanaka pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi scheme. The first non-dimen ¼ xh qm =Em according to various meshes sional frequencies x are listed in Table 5 and then compared with that from the exact solution by Vel and Batra [1], the solutions based on HSDT proposed by Ferreira et al. [15] using the collocation multiquadric radial basis functions and Qian et al. [29] using meshless local Petrov–Galerkin (MLPG) method. Considering the coupling of the

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

12

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

Table 6 The first ten non-dimensional frequencies of simply supported Al/ZrO2 FG plate with a/h = 5. n

Method or authors

Mode number Mode1

Mode2

Mode3

Mode4

Mode 5

Mode 6

Mode 7

Mode 8

Mode 9

Mode 10

0

IGA Meshless method [15] MLPG [29]

0.2461 0.2457 0.2469

0.4539 0.4483 0.4535

0.4539 0.4484 0.4535

0.5385 0.5395 0.5441

0.5385 0.5395 0.5441

0.6419 0.647 0.6418

0.7791 0.7809 0.7881

0.9078 0.899 0.9076

0.9078 0.8994 0.9326

0.9210 0.9209 0.9354

1

IGA Meshless method [15] MLPG [29]

0.2185 0.2188 0.2152

0.4118 0.399 0.4114

0.4118 0.3992 0.4114

0.4794 0.4779 0.4761

0.4794 0.4779 0.4761

0.5823 0.5759 0.582

0.6948 0.6897 0.6914

0.8219 0.8001 0.8192

0.8219 0.8163 0.8217

0.8232 0.8118 0.8242

2

IGA Meshless method [15] MLPG [29]

0.2190 0.2188 0.2153

0.4039 0.3991 0.4034

0.4039 0.3992 0.4034

0.4768 0.4779 0.472

0.4768 0.4779 0.472

0.5711 0.5759 0.5709

0.6878 0.6897 0.6817

0.8075 0.8001 0.8056

0.8075 0.8004 0.8105

0.8119 0.8118 0.9022

5

IGA Meshless method [15] MLPG [29]

0.2216 0.2215 0.2194

0.3967 0.3921 0.3964

0.3967 0.3922 0.3964

0.4783 0.4794 0.476

0.4783 0.4795 0.476

0.5611 0.5658 0.5611

0.6865 0.6884 0.6832

0.7934 0.7863 0.7928

0.7934 0.7866 0.8053

0.8084 0.8081 0.8099

Metal

IGA Meshless method [15] MLPG [29]

0.2113 0.2111 0.2122

0.3897 0.3852 0.3897

0.3897 0.3853 0.3897

0.4623 0.4636 0.4675

0.4623 0.4636 0.4675

0.5511 0.556 0.5517

0.6688 0.671 0.6772

0.7793 0.7725 0.7615

0.7793 0.7728 0.7799

0.7906 0.7913 0.8013

Table 7 pffiffiffiffiffiffiffiffiffiffiffiffi  ¼ 10xh qc =Ec of simply supported square Al/Al2O3 FG The first natural frequency x plate. a/h

5

10

20

Theory

Method

n 0

0.5

1

4

10

FSDT

Exact NS-DSG3 ES-DSG3 DSG3 MITC4 IGA

2.112 2.114 2.127 2.133 2.118 2.1398

1.805 1.805 1.816 1.822 1.808 1.8259

1.631 1.630 1.640 1.644 1.632 1.6486

1.397 1.394 1.403 1.407 1.397 1.4135

1.324 1.323 1.331 1.334 1.325 1.3425

HSDT

2D-HDT IGA

2.121 2.1125

1.819 1.8068

1.640 1.6306

1.383 1.3773

1.306 1.2995

FSDT

Exact NS-DSG3 ES-DSG3 DSG3 MITC4 IGA

0.577 0.578 0.582 0.583 0.579 0.5794

0.490 0.491 0.494 0.495 0.491 0.4917

0.442 0.443 0.445 0.447 0.443 0.4433

0.382 0.383 0.385 0.386 0.383 0.3837

0.366 0.366 0.368 0.369 0.366 0.3673

HSDT

2D-HDT IGA

0.578 0.5769

0.492 0.4900

0.443 0.4418

0.381 0.3804

0.364 0.3634

FSDT

Exact NS-DSG3 ES-DSG3 DSG3 MITC4 IGA

0.148 0.148 0.149 0.150 0.148 0.1482

0.125 0.126 0.127 0.127 0.126 0.1255

0.113 0.113 0.114 0.114 0.113 0.1131

0.098 0.098 0.099 0.099 0.098 0.0982

0.094 0.094 0.095 0.095 0.094 0.0943

HSDT

IGA

0.1480

0.1254

0.1130

0.0980

0.0940

high-order terms with negative c given in Eq. (5) the HSDT model gets the softer stiffness matrix than FSDT one, thus result of the former is also lower than that of the latter. As compared to exact values [1] the present method gains the lowest error which is less than 0.4%. In addition, the IGA gives the supper convergence in dynamics analysis with closely zero discrepancy between meshing 5  5 and 15  15. Table 6 reveals the first ten non-dimensional frequencies with the meshing of 11  11 NURBS elements of the moderate Al/ZrO2  1 plate. The present computed values agree well with the literature [15,29] via various power indexes. Next, the rule of mixture is applied in the investigation of free vibration of the Al/Al2O3 plate under the change of length to thickness ratios. Table 7 shows the comparison between present results and those from HSDT by [37] and FSDT [11,12,38,39]. Some remarks are given here: (1) with the same used basic function, because of stiffer stiffness matrix, FSDT gains the slightly lar-

ger frequency parameter than HSDT; (2) the thinner the plate becomes, the lower discrepancy between FSDT and HSDT’s solution is; (3) As n 6 1 among the numerical methods, the present method gains the closest value to analytical solution [39], while as n > 1, the present results agree well with that from 2D higher-order deformation theory (2D-HDT) approach by Matsunaga [37]; (4) The frequency parameter reduces according to the decrease of length to thick a/h ratio because of weakness of the plate.

4.2.2. Circular plate pffiffiffiffiffiffiffiffiffiffiffiffi  ¼ 100xh qc =Ec of In this section, the normalized frequency x the clamped circular Al/Al2O3 plate with radius R = 0.5 are tested. The effective Young’s modulus Eeff and density qeff of FG plate are calculated according to Eq. (36). Herein, the effects of the power index and radius to thickness ratio R/h are investigated. The obtained results as shown in Table 8 agree well with those by three methods: semi-analytical solution [40], a finite element method from ABAQUS software package and a method named the uncoupled model (UM), formerly proposed by Ebrahimi et al. [41]. To close this section, the first six mode shapes of circular plate are plotted in Fig. 16.

4.3. Buckling analysis Let’s consider a clamped plate of radius R and thickness h subjected to a uniform radial pressure p0 as shown in Fig. 17. The dependent of buckling load parameters of Al/Al2O3 plates on radius to thickness ratio are revealed in Fig. 18. It is concluded that pcr increases with the decrease of radius to thickness ratio. And the homogeneous ceramic plate (n = 0) with highest Young’s modulus value gets the highest buckling load parameter. Those obtained results are then compared to closed form solution [43] summarized in Table 9. It is seen that the present results are slightly higher than that by M.M. Najafizadeh using the same HSDT model. The discrepancy between them reduces from 0.27% to 0.16% as R/h ranges from 0.01 to 0.1. The last, buckling response of Al/ZrO2  2 thick plate which material properties are computed following to Eq. (36) are investigated. Table 10 shows that the present results are well compared to those analytical solution in the published literatures [42,44] using HSDT assumption. Fig. 19 illustrates the variation of the critical load versus the power index. It is observed that the critical

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

13

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx Table 8 pffiffiffiffiffiffiffiffiffiffiffiffi  ¼ 100xh qc =Ec of circular Al/Al2O3 FG plate with clamped edge. The first eight frequencies x h/R

Method

Mode number Mode 2

Mode 3

Mode 4

Mode 5

Mode 6

Mode 7

Mode 8

0.01

Ref FEM UM IGA

0.0236 0.0234 0.0257 0.0236

0.0491 0.0486 0.0535 0.0492

0.0805 0.0798 0.0877 0.0807

0.0918 0.0909 0.1000 0.0924

0.1178 0.1167 0.1283 0.1191

0.1404 0.1391 0.1529 0.1431

0.1607 0.1592 0.1751 0.1643

0.1951 0.1933 0.2126 0.1991

0.1

Ref FEM UM IGA

2.3053 2.2888 2.5038 2.3076

4.6934 4.6661 5.0831 4.7005

7.5146 7.4808 8.1156 7.5318

8.5181 8.4829 9.1931 8.5380

10.7128 10.6776 11.5376 10.7483

12.6197 12.5877 13.5743 12.6636

14.2324 14.2025 15.2879 14.2925

16.9838 16.9583 18.2114 17.0520

0.2

Ref FEM UM IGA

8.6535 8.6403 9.3162 8.6787

16.7666 16.7890 17.9164 16.8595

25.6486 25.7661 27.2480 25.8479

28.7574 28.9152 30.4998 29.0092

34.0756 34.1893 – 34.0581

35.0981 35.3618 37.1197 35.4875

39.4394 39.4169 – 39.4177

40.5889 40.9538 42.8001 41.0759

Mode 1

Fig. 16. The six first mode shapes of clamped circular Al/Al2O3 plate with h/R = 0.1.

Fig. 17. The circular plate subjected to uniform radial compression.

cr ¼ pcr R2 =ðEm h3 Þ of the Al/Al2O3 plate via Fig. 18. The normalized buckling load p radius to thickness ratio R/h.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

14

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx Table 9 The critical buckling load pcr = p0R2/(Emh3) of thin and moderate circular Al/Al2O3 plates via h/R ratio and power index n. Power index

Method

h/R 0.01

0.02

0.05

0.08

0.1

0

HSDT FSDT CPT IGA (p = 3)

7.2930 7.2958 7.2977 7.3119

7.2838 7.2899 7.2977 7.2958

7.2204 7.2490 7.2977 7.2321

7.1054 7.1743 7.2977 7.1168

7.0025 7.1067 7.2977 7.0137

0.5

HSDT FSDT CPT IGA (p = 3)

4.7279 4.7296 4.7307 4.7406

4.7227 4.7262 4.7307 4.7303

4.6865 4.7026 4.7307 4.6940

4.6209 4.6593 4.7307 4.6282

4.5619 4.6201 4.7307 4.5691

2

HSDT FSDT CPT IGA (p = 3)

2.8367 2.8377 2.8384 2.8444

2.8334 2.8358 2.8384 2.8383

2.8111 2.8222 2.8384 2.8159

2.7705 2.7974 0.2838 2.7752

2.7341 2.7748 2.8384 2.7387

5

HSDT FSDT CPT IGA (p = 3)

2.3985 2.3997 2.4004 2.4046

2.3948 2.3977 2.4004 2.3991

2.3691 2.3839 2.4004 2.3733

2.3229 2.3587 0.2400 2.3269

2.2817 2.3358 2.4004 2.2857

Table 10 cr ¼ pcr R2 =Dm of clamped thick circular Al/ZrO2  2 plate. Comparison of the buckling load parameter p Power index

Method

h/R 0.1

0.2

0.25

0.3

0

TST UTST IGA (p = 3)

14.089 14.089 14.1089

12.574 12.575 12.5914

11.638 11.639 11.6540

10.67 10.67 10.6842

0.5

TST UTST IGA (p = 3)

19.411 19.413 19.4391

17.311 17.31 17.3327

16.013 16.012 16.0334

14.672 14.672 14.6910

2

TST UTST IGA (p = 3)

23.074 23.075 23.1062

20.803 20.805 20.8319

19.377 19.378 19.4033

17.882 17.881 17.9060

5

TST UTST IGA (p = 3)

25.439 25.442 25.4743

22.971 22.969 22.9992

21.414 21.412 21.4407

19.78 19.778 19.8043

10

TST UTST IGA (p = 3)

27.133 27.131 27.1684

24.423 24.422 24.4542

22.725 22.725 22.7536

20.948 20.949 20.9750

5. Conclusions In this paper, a novel plate formulation using NURBS-based isogeometric approach and HSDT has been applied for static, dynamic and buckling analysis of functionally graded materials plates. The effective material properties at arbitrary point in the FGM plates are computed either the rule of mixture or the Mori–Tanaka homogenization technique. The present method allows us to achieve easily the smoothness with arbitrary continuity order compared with the traditional FEM. As a result, the present NURBS formulation naturally fulfills the C1-continuity of HSDT plates. The results of the present method are in very excellent agreement with several other published solutions in the literature. Furthermore, some additional discussions can be drawn from above numerical examples: cr ¼ pcr R2 =Dm of the Al/ZrO2  2 plate via Fig. 19. The normalized buckling load p power index n.

cr ¼ pcr R2 =Dm changes rapidly as n ranges buckling load parameter p from 0 to 10. But, as n > 10 it is slightly independent on the increment of power index.

– Using isogeometric approach, not only the plate geometry is described exactly (e.g. circular plate), but also the obtain results match very well with published others. Moreover, the exact geometry representation is inherited from the first coarsest mesh and preserved during the refinement process without the need for additional communication with the CAD program.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

– The Mori–Tanaka homogeneous scheme considering the interactions among the constituents gets the lower property feature than the rule of mixtures. Hence, the former gains larger displacements, it is rather suitable in design FGM plates for safety. – HSDT model considering the coupling of the high-order terms with negative value c given in Eq. (5) gains the softer stiffness matrix than FSDT one. Thus, in static problems, the former gets higher results which are closer to 3D deformation solutions than the latter. In addition, with third shear deformation assumption, the present method gains the shear stress which is in form of parabolically curved distribution and satisfies the traction free conditions at top and bottom surfaces of the plates. In free vibration and mechanical buckling problems, discrepancies between HSDT and FSDT’s results are very small. In author’s opinion, to simply calculate, we can use FSDT model for those problems of FGM plates. – Axial stress through FGM plate thickness is distributed in form taking advantage of material characteristic with increase of compression at the top (ceramic surface) and decrease of tension at the bottom (metal surface), respectively. In this paper, only the third order shear deformation theory is utilized. However, the IGA can be extended for other theories such as: 3D model [45], a quasi-3D sinusoidal shear deformation theory [46], higher-order shear deformation theory considering stretching effect [47] and some plate models with trigonometric approximation through the thickness coordinate [48–50]. It is believed that the proposed method can be very promising to provide the effectively computational tool for practical problems. For purpose of comparison with the existing solution formulations, we only restrict to our numerical study to the benchmarks of simple geometries such as square and circle plates. However, the IGA based on NURBS is straightforward for describing more complex domains with many holes and chamfers by multi-patch technique [25] or finite cell method [51]. This problem, which is very useful for analysis of real structures, will be studied in future work. Acknowledgement This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.02-2012.17. The support is gratefully acknowledged. References [1] Vel SS, Batra RC. Exact solution for thermoelastic deformations of functionally graded thick rectangular plates. AIAA J 2002;40:1021–33. [2] Koizumi M. FGM activities in Japan. Compos Part B: Eng 1997;28(1–2):1–4. [3] Fukui Y, Yamanaka N. Elastic analysis for thick-walled tubes of functionally graded material subjected to internal pressure. Int J Jpn Soc Mech Eng A 1992;35:379–85. [4] Obata Y, Noda N. Transient thermal stresses in a plate of functionally gradient material. Ceram Trans Funct Grad Mater 1993;34:403–10. [5] Obata Y, Noda N. Optimum material design for functionally gradient material plate. Arch Appl Mech 1996;66:581–9. [6] Reddy JN. Analysis of functionally graded plates. Int J Numer Methods Eng 2000;47:663–84. [7] Frohlich H, Sack R. Theory of the rheological properties of dispersions. Proc R Soc Lond A Math Phys Sci 1946;185:415–30. [8] Hill R. A self-consistent mechanics of composite materials. J Mech Phys Solids 1965;13:213–22. [9] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;21:571–4. [10] Jiang B, Batra RC. Effective properties of a piezocomposite containing shape memory alloy and inert inclusions. Continuum Mech Thermodynam 2002;14:87–111. [11] Nguyen-Xuan H, Tran Loc V, Nguyen-Thoi T, Vu-Do HC. Analysis of functionally graded plates using an edge-based smoothed finite element method. Compos Struct 2011;93:3019–39. [12] Nguyen-Xuan H, Tran Loc V, Thai Chien H, Nguyen-Thoi T. Analysis of functionally graded plates by an efficient finite element method with nodebased strain smoothing. Thin-Wall Struct 2012;54:1–18.

15

[13] Prakash T, Ganapathi M. Asymmetric flexural vibration and thermo elastic stability of FGM circular plates using finite element method. Composites: Part B 2006;37:642–9. [14] Ferreira AJM, Batra RC, Roque CMC, Qian LF, Martins PALS. Static analysis of functionally graded plates using third-order shear deformation theory and a meshless method. Compos Struct 2005;69:449–57. [15] Ferreira AJM, Batra RC, Roque CMC, Qian LF, Jorge RMN. Natural frequencies of functionally graded plates by a meshless method. Compos Struct 2006;75:593–600. [16] Thai Chien H, Tran Loc V, Tran Dung T, Nguyen-Thoi T, Nguyen-Xuan H. Analysis of laminated composite plates using higher-order shear deformation plate theory and node-based smoothed discrete shear gap method. Appl Math Modell 2012;36:5657–77. [17] Tran Loc V, Nguyen-Thoi T, Thai Chien H, Nguyen-Xuan H. An edge-based smoothed discrete shear gap method (ES-DSG) using the C0-type higher-order shear deformation theory for analysis of laminated composite plates. Mech Adv Mater Struct, 2012, in press. [18] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194(39–41):4135–95. [19] Cottrell JA, Hughes TJR, Bazilevs Y. Isogeometric analysis, towards integration of CAD and FEA. Wiley; 2009. [20] lguedj TE, Bazilevs Y, Calo V, Hughes T. B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higherorder NURBS elements. Comput Methods Appl Mech Eng 2008;197:2732–62. [21] Cottrell JA, Reali A, Bazilevs Y, Hughes TJR. Isogeometric analysis of structural vibrations. Comput Methods Appl Mech Eng 2006;195(41–43):5257–96. [22] Beirão da Veiga L, Buffa A, Lovadina C, Martinelli M, Sangalli G. An isogeometric method for the Reissner–Mindlin plate bending problem. Comput Methods Appl Mech Eng 2012;209–212:45–53. [23] Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. Isogeometric shell analysis: The Reissner–Mindlin shell. Comput Methods Appl Mech Eng 2010;199(5– 8):276–89. [24] Kiendl J, Bletzinger KU, Linhard J, Wchner R. Isogeometric shell analysis with Kirchhoff–Love elements. Comput Methods Appl Mech Eng 2009;198(49– 52):3902–14. [25] Kiendl J, Bazilevs Y, Hsu MC, Wchner R, Bletzinger KU. The bending strip method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches. Comput Methods Appl Mech Eng 2010;199(37– 40):2403–16. [26] Wall WA, Frenzel MA, Cyron C. Isogeometric structural shape optimization. Comput Methods Appl Mech Eng 2008;197(33–40):2976–88. [27] Nguyen-Thanh N, Nguyen-Xuan H, Bordas S, Rabczuk T. Isogeometric analysis using polynomial splines over hierarchical for two-dimensional elastic solids. Comput Methods Appl Mech Eng 2011;200:1892–908. [28] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wüchner R, Bletzinger KU, Bazilevs Y, et al. Rotation free isogeometric thin shell analysis using PHTsplines. Comput Methods Appl Mech Eng 2011;200(47–48):3410–24. [29] Qian LF, Batra RC, Chen LM. Static and dynamic deformations of thick functionally graded elastic plate by using higher-order shear and normal deformable plate theory and meshless local Petrov–Galerkin method. Compos Part B: Eng 2004;35(6–8):685–97. [30] Reddy JN. A simple higher-order theory for laminated composite plates. J Appl Mech 1984;51:745–52. [31] Reddy JN. Mechanics of laminated composite plates – theory and analysis. New York: CRC Press; 1997. [32] Sankara CA, Igengar NGR. A C0 element for free vibration analysis of laminated composite plates. J Sound Vib 1996;19:721–38. [33] Thai Chien H, Nguyen-Xuan H, Nguyen-Thanh N, Le T-H, Nguyen-Thoi T, Rabczuk T. Static, free vibration, and buckling analysis of laminated composite Reissner–Mindlin plates using NURBS-based isogeometric approach. Int J Numer Methods Eng 2012. http://dx.doi.org/10.1002/nme.4282. [34] Gilhooley DF, Batra RC, Xiao JR, McCarthy MA, Gillespie JW. Analysis of thick functionally graded plates by using higher-order shear and normal deformable plate theory and MLPG method with radial basis functions. Comput Struct 2007;80:539–52. [35] Li XY, Ding HJ, Chen WQ. Elasticity solutions for a transversely isotropic functionally graded circular plate subject to an axisymmetric transverse load qrk. Int J Solids Struct 2008;45:191–210. [36] Reddy JN, Wang CM, Kitipornchai S. Axisymmetric bending of functionally graded circular and annular plates. Eur J Mech A Solids 1999;18:185–99. [37] Matsunaga H. Free vibration and stability of functionally graded plates according to a 2-D higher-order deformation theory. Compos Struct 2008;82:499–512. [38] Zhao X, Lee YY, Liew KM. Free vibration analysis of functionally graded plates using the element-free kp-Ritz method. J Sound Vib 2003;319:918–39. [39] Hosseini-Hashemi Sh, Fadaee M, Atashipour SR. A new exact analytical approach for free vibration of Reissner–Mindlin functionally graded rectangular plates. Int J Mech Sci 2011;53:11–22. [40] Hosseini-Hashemi Sh, Fadaee M, Es’haghi M. A novel approach for in-plane/ out-of-plane frequency analysis of functionally graded circular/annular plates. Int J Mech Sci 2010;52:1025–35. [41] Ebrahimi F, Rastgoo A, Atai AA. A theoretical analysis of smart moderately thick shear deformable annular functionally graded plate. Eur J Mech A Solids 2009;28:962–73.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045

16

L.V. Tran et al. / Composites: Part B xxx (2013) xxx–xxx

[42] Saidi AR, Rasouli A, Sahraee S. Axisymmetric bending and buckling analysis of thick functionally graded circular plates using unconstrained third-order shear deformation plate theory. Compos Struct 2009;89:110–9. [43] Najafizadeh MM, Heydari HR. An exact solution for buckling of functionally graded circular plates based on higher order shear deformation plate theory under uniform radial compression. Int J Mech Sci 2008;50:603–12. [44] Ma LS, Wang TJ. Relationship between axisymmetric bending and buckling solutions of FGM circular plates based on third-order plate theory and classical plate theory. Int J Solids Struct 2004;41:85–101. [45] Vaghefi R, Baradaran GH, Koohkan H. Three-dimensional static analysis of thick functionally graded plates by using meshless local Petrov–Galerkin (MLPG) method. Eng Anal Bound Elem 2010;34:564–73. [46] Neves AMA, Ferreira AJM, Carrera E, Roque CMC, Cinefra M, Jorge RMN, et al. A quasi-3D sinusoidal sheared formation theory for the static and free vibration analysis of functionally graded plates. Composites: Part B 2012;43:711–25. [47] Mantari JL, Guedes Soares C. A novel higher-order shear deformation theory with stretching effect for functionally graded plates. Compos Part B: Eng 2013;45:268–81.

[48] Mantari JL, Oktem AS, Guedes Soares C. Bending response of functionally graded plates by using a new higher order shear deformation theory. Compos Struct 2012;94:714–23. [49] Zenkour AM. The refined sinusoidal theory for FGM plates on elastic foundations. Int J Mech Sci 2009;51(11–12):869–80. [50] Ait Atmane H, Tounsi A, Mechab I, Adda Bedia EA. Free vibration analysis of functionally graded plates resting on Winkler–Pasternak elastic foundations using a new shear deformation theory. Int J Mech Mater Des 2010;6(2):113–21. [51] Schillinger D, Ruess M, Zander N, Bazilevs Y, Düster A, Rank E. Small and large deformation analysis with the p- and B-spline versions of the Finite Cell Method. Comput Mech 2012;50:445–78. [52] Thai Chien H, Nguyen-Xuan H, Bordas SPA, Nguyen-Thanh N, Rabczuk T. Isogeometric analysis of laminated composite plates using the higher-order shear deformation theory. Mech Adv Mater Struct, 2012, in press.

Please cite this article in press as: Tran LV et al. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites: Part B (2013), http://dx.doi.org/10.1016/j.compositesb.2013.02.045