# Analysis of Functionally Graded Material Plates Using Triangular Elements with Cell-Based Smoothed Discrete Shear Gap Method

Authors: S. Natarajan; A. J. M. Ferreira; S. Bordas; E. Carrera; M. Cinefra; A. M. Zenkour
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A cell-based smoothed finite element method with discrete shear gap technique is employed to study the static bending, free vibration, and mechanical and thermal buckling behaviour of functionally graded material (FGM) plates. The plate kinematics is based on the first-order shear deformation theory and the shear locking is suppressed by the discrete shear gap method. The shear correction factors are evaluated by employing the energy equivalence principle. The material property is assumed to be temperature dependent and graded only in the thickness direction. The effective properties are computed by using the Mori-Tanaka homogenization method. The accuracy of the present formulation is validated against available solutions. A systematic parametric study is carried out to examine the influence of the gradient index, the plate aspect ratio, skewness of the plate, and the boundary conditions on the global response of the FGM plates. The effect of a centrally located circular cutout on the global response is also studied.


---

## Body

## 1. Introduction

With the rapid advancement of engineering, there is an increasing demand for new materials which suit the harsh working environment without losing their mechanical, thermal, or electrical properties. Engineered materials such as the composite materials are used due to their excellent strength-to-weight and stiffness-to-weight ratios and their possibility of tailoring the properties in optimizing their structural response. But due to the abrupt change in material properties from matrix to fibre and between the layers, these materials suffer from premature failure or from the decay in the stiffness characteristics because of delaminations and chemically unstable matrix and lamina adhesives. On the contrary, another class of materials, called the functionally graded materials (FGMs), are made up of mixture of ceramics and metals and are characterized bysmooth and continuous transition in properties from one surface to another [1]. As a result, FGMs are preferred over the laminated composites for structural integrity. The FGMs combine the best properties of the ceramics and the metals and this has attracted the researchers to study the characteristics of such structures.Background. The tunable thermomechanical property of the FGM has attracted researchers to study the static and the dynamic behaviour of structures made of FGM under mechanical [2–5] and thermal loading [6–12]. Praveen and Reddy [7] and Reddy and Chin [13] studied the thermoelastic response of ceramic-metal plates using first-order shear deformation theory coupled with the 3D heat conduction equation. Their study concluded that the structures made up of FGM with ceramic rich side exposed to elevated temperatures are susceptible to buckling due to the through-thickness temperature variation. The buckling of skewed FGM plates under mechanical and thermal loads was studied [9, 14] by employing the first-order shear deformation theory and by using the shear flexible quadrilateral element. Efforts have also been made to study the mechanical behaviour of FGM plates with geometrical imperfection [15]. Saji et al. [16] have studied thermal buckling of FGM plates with material properties dependent on both the composition and temperature. They found that the critical buckling temperatures decreased when material properties are considered to be a function of temperature as compared to the results obtained where the material properties are assumed to be independent of temperature. Ganapathi and Prakash [9] have studied the buckling of FGM skewed plate under thermal loading. Efforts have also been made to study the mechanical behaviour of FGM plates with geometrical imperfection [15]. More recently, refined models have been adopted to study the characteristics of FGM structures [17–19].Existing approaches in the literature to study plate and shell structures made up of FGMs use finite element method (FEM) based on Lagrange basis functions [14] and meshfree methods [20, 21], and recently Valizadeh et al. [22] used nonuniform rational B-splines based FEM to study the static and dynamic characteristics of FGM plates in thermal environment. Loc et al. [23] employed isogeometric finite element method to study thermal buckling of functionally graded plates. Even with these different approaches, the plate elements suffer from shear locking phenomenon and different techniques were proposed to alleviate the shear locking phenomenon. Another set of methods have emerged to address the shear locking in the FEM. By incorporating the strain smoothing technique into the finite element method (FEM), Liu et al. [24] have formulated a series of smoothed finite element methods (SFEM), named as cell-based SFEM (CS-FEM) [25, 26], node-based SFEM [27], edge-based SFEM [28], face-based SFEM [29], and α-FEM [30]. And recently, edge based imbricate finite element method (EI-FEM) was proposed in [31] that shares common features with the ES-FEM. As the SFEM can be recast within a Hellinger-Reissner variational principle, suitable choices of the assumed strain/gradient space provide stable solutions. Depending on the number and geometry of the subcells used, a spectrum of methods exhibiting a spectrum of properties is obtained. Interested readers are referred to the literature [24, 25] and references therein. Nguyen-Xuan et al. [32] employed CS-FEM for Mindlin-Reissner plates. The curvature at each point is obtained by a nonlocal approximation via a smoothing function. From the numerical studies presented, it was concluded that the CS-FEM technique is robust, computationally inexpensive, free of locking, and most importantly insensitive to mesh distortions. The SFEM was extended to various problems such as shells [33], heat transfer [34], fracture mechanics [35], and structural acoustics [36] among others. In [37], CS-FEM has been combined with the extended FEM to address problems involving discontinuities. The above list is no way comprehensive, and interested readers are referred to the literature and references therein and to a recent review paper by Jha et al. [38] on FGM plates.Approach. In this paper, we study the static and the dynamic characteristics of FGM plates by using a cell-based smoothed finite element method with discrete shear gap technique [39]. Three-noded triangular element is employed in this study. The effect of different parameters, namely, the material gradient index, the plate aspect ratio, the plate slenderness ratio, and the boundary condition on the global response of FGM plates, is numerically studied. The effect of centrally located circular cutout is also studied. The present work focuses on the computational aspects of the governing equations and hence the attention has been restricted to Mindlin-Reissner plate theory. It is noted that the extension to higher order theories is possible.Outline. The paper is organized as follows: the next section will give an introduction to FGM and a brief overview of Mindlin-Reissner plate theory. Section 3 presents an overview of the cell-based smoothed finite element method with discrete shear gap technique. The efficiency of the present formulation, numerical results, and parametric studies are presented in Section 4, followed by concluding remarks in the last section.


## 2. Theoretical Background

### 2.1. Mindlin-Reissner Plate Theory

The Mindlin-Reissner plate theory, also known as the first-order shear deformation theory (FSDT), takes into account the shear deformation through the thickness. Using the Mindlin formulation, the displacementsu, v, and w at a point (
x
,
y
,
z
) in the plate (see Figure 1) from the medium surface are expressed as functions of the midplane displacements u
o, v
o, and w
o and independent rotations θ
x and θ
y of the normal in y
z and x
z planes, respectively, as follows:
(1)
u
(
x
,
y
,
z
,
t
)
=
u
o
(
x
,
y
,
t
)
+
z
θ
x
(
x
,
y
,
t
)
,
v
(
x
,
y
,
z
,
t
)
=
v
o
(
x
,
y
,
t
)
+
z
θ
y
(
x
,
y
,
t
)
,
w
(
x
,
y
,
z
,
t
)
=
w
o
(
x
,
y
,
t
)
,
where t is the time. The strains in terms of midplane deformation can be written as follows:
(2)
ε
=
{
ε
p
0
}
+
{
z
ε
b
ε
s
}
.
The midplane strains ε
p, the bending strain ε
b, and the shear strain ε
s in (2) are written as follows:
(3)
ε
p
=
{
u
o
,
x
v
o
,
y
u
o
,
y
+
v
o
,
x
}
,
ε
b
=
{
θ
x
,
x
θ
y
,
y
θ
x
,
y
+
θ
y
,
x
}
,
ε
s
=
{
θ
x
+
w
o
,
x
θ
y
+
w
o
,
y
}
,
where the subscript “comma” represents the partial derivative with respect to the spatial coordinate succeeding it. The membrane stress resultants N and the bending stress resultants M can be related to the membrane strains ε
p and bending strains ε
b through the following constitutive relations:
(4)
N
=
{
N
x
x
N
y
y
N
x
y
}
=
A
ε
p
+
B
ε
b
-
N
th
,
M
=
{
M
x
x
M
y
y
M
x
y
}
=
B
ε
p
+
D
b
ε
b
-
M
th
,
where the matrices A
=
A
i
j
,
B
=
B
i
j, and D
b
=
D
i
j
;
(
i
,
j
=
1,2
,
6
) are the extensional and bending-extensional coupling and bending stiffness coefficients and are defined as follows:
(5)
{
A
i
j
,
B
i
j
,
D
i
j
}
=
∫
-
h
/
2
h
/
2
‍
Q
¯
i
j
{
1
,
z
,
z
2
}
d
z
.
Similarly, the transverse shear force Q
=
{
Q
x
z
,
Q
y
z
} is related to the transverse shear strains ε
s through the following equation:
(6)
Q
i
j
=
E
i
j
ε
s
,
where E
=
E
i
j
=
∫
-
h
/
2
h
/
2
‍
Q
¯
i
j
υ
i
υ
j
d
z
;
(
i
,
j
=
4,5
) are the transverse shear stiffness coefficients, and υ
i, υ
j is the transverse shear coefficient for nonuniform shear strain distribution through the plate thickness. The stiffness coefficients Q
¯
i
j are defined as
(7)
Q
¯
11

=

Q
¯
22
=
E
(
z
)
1
-
ν
2
;
Q
¯
12
=
ν
E
(
z
)
1
-
ν
2
;
Q
¯
16
=
Q
¯
26
=
0
,
Q
¯
44
=
Q
¯
55
=
Q
¯
66
=
E
(
z
)
2
(
1
+
ν
)
,
where the modulus of elasticity E
(
z
) and Poisson’s ratio ν are given by (20). The thermal stress resultant N
th and the moment resultant M
th are
(8)
N
th
=
{
N
x
x
th
N
y
y
th
N
x
y
th
}
=
∫
-
h
/
2
h
/
2
Q
¯
i
j
α
(
z
,
T
)
{
1
1
0
}
△
T
(
z
)
d
z
,
M
th
=
{
M
x
x
th
M
y
y
th
M
x
y
th
}
=
∫
-
h
/
2
h
/
2
Q
¯
i
j
α
(
z
,
T
)
{
1
1
0
}
△
T
(
z
)
z
d
z
,
where the thermal coefficient of expansion α
(
z
,
T
) is given by (21), △
T
(
z
)
=
T
(
z
)
-
T
o is the temperature rise from the reference temperature, and T
o is the temperature at which there are no thermal strains. The strain energy function U is given by
(9)
U
(
δ
)
=
1
2
∫
Ω
‍
{
ε
p
T
A
ε
p
+
ε
p
T
B
ε
b
+
ε
b
T
B
ε
p
k
k
k
k
k
k
l
l
+
ε
b
T
D
ε
b
+
ε
s
T
E
ε
s
-
ε
b
T
N
th
-
ε
b
T
M
th
}
d
Ω
,
where δ
=
{
u
,
v
,
w
,
θ
x
,
θ
y
} is the vector of the degree of freedom associated with the displacement field in a finite element discretization. Following the procedure given in [40], the strain energy function U given in (9) can be rewritten as
(10)
U
(
δ
)
=
1
2
δ
T
K
δ
,
where K is the linear stiffness matrix. The kinetic energy of the plate is given by
(11)
T
(
δ
)
=
1
2
∫
Ω
‍
{
p
(
u
˙
o
2
+
v
˙
o
2
+
w
˙
o
2
)
+
I
(
θ
˙
x
2
+
θ
˙
y
2
)
}
d
Ω
,
where =
∫
-
h
/
2
h
/
2
‍
ρ
(
z
)
d
z
,
I
=
∫
-
h
/
2
h
/
2
‍
z
2
ρ
(
z
)
d
z, and ρ
(
z
) is the mass density that varies through the thickness of the plate. When the plate is subjected to a temperature field, this in turn results in in-plane stress resultants, N
th. The external work due to the in-plane stress resultants developed in the plate under a thermal load is given by
(12)
V
(
δ
)
=
∫
Ω
{
1
2
[
N
x
x
th
w
,
x
2
+
N
y
y
th
w
,
y
2
+
2
N
x
y
th
w
,
x
w
,
y
]
kkkl
+
h
2
24
[
N
x
x
th
(
θ
x
,
x
2
+
θ
y
,
x
2
)
+
N
y
y
2
(
θ
x
,
y
2
+
θ
y
,
y
2
)
kkkl
hh
+
2
N
x
y
th
(
θ
x
,
x
θ
x
,
y
+
θ
y
,
x
θ
y
,
y
)
]
1
2
}
d
Ω
.
Substituting (9)–(12) in Lagrange’s equation of motion, one obtains the following finite element equations:Static bending:(13)
K
δ
=
F
.
Free vibration:(14)
M
δ
¨
+
(
K
+
K
G
)
δ
=
0
.
Buckling analysis:
Mechanical buckling (prebuckling deformations are assumed to be zero or negligible in the analysis (including those coming from in-plane and out-of-plane coupling related to FGM and temperature variation through the thickness of the plate)):(15)
(
K
+
λ
M
K
G
)
δ
=
0
.
Thermal buckling:(16)
(
K
+
λ
T
K
G
)
δ
=
0
, where λ
M includes the critical value applied to in-plane mechanical loading, λ
T is the critical temperature difference, and K and K
G are the linear stiffness and geometric stiffness matrices, respectively. The critical temperature difference is computed using a standard eigenvalue algorithm.(a) Coordinate system of a rectangular FGM plate and (b) coordinate system of a skew plate.
(a)
(b)


### 2.2. Functionally Graded Material

A rectangular plate made of a mixture of ceramic and metal is considered with the coordinatesx
,
y along the in-plane directions and z along the thickness direction (see Figure 1). The material on the top surface (
z
=
h
/
2
) of the plate is ceramic rich and is graded to metal at the bottom surface of the plate (
z
=
-
h
/
2
) by a power law distribution. The effective properties of the FGM plate can be computed by using the rule of mixtures or by employing the Mori-Tanaka homogenization scheme. Let V
i
(
i
=
c
,
m
) be the volume fraction of the phase material. The subscripts c and m refer to ceramic and metal phases, respectively. The volume fraction of ceramic and metal phases is related by V
c
+
V
m
=
1, and V
c is expressed as
(17)
V
c
(
z
)
=
(
2
z
+
h
2
h
)
n
,
where n is the volume fraction exponent (
n
≥
0
), also known as the gradient index. The variation of the composition of ceramic and metal is linear for n
=
1, the value of n
=
0 represents a fully ceramic plate, and any other value of n yields a composite material with a smooth transition from ceramic to metal.Mori-Tanaka Homogenization Method. Based on the Mori-Tanaka homogenization method, the effective Young’s modulus and Poisson’s ratio are computed from the effective bulk modulus K and the effective shear modulus G as [41]
(18)
K
eff
-
K
m
K
c
-
K
m
=
V
c
1
+
V
m
(
3
(
K
c
-
K
m
)
/
(
3
K
m
+
4
G
m
)
)
,
G
eff
-
G
m
G
c
-
G
m
=
V
c
1
+
V
m
(
(
G
c
-
G
m
)
/
(
G
m
+
f
1
)
)
,
where
(19)
f
1
=
G
m
(
9
K
m
+
8
G
m
)
6
(
K
m
+
2
G
m
)
.
The effective Young’s modulus E
eff and Poisson’s ratio ν
eff can be computed from the following relations:
(20)
E
eff
=
9
K
eff
G
eff
3
K
eff
+
G
eff
,
ν
eff
=
3
K
eff
-
2
G
eff
2
(
3
K
eff
+
G
eff
)
.
The effective mass density ρ
eff is computed using the rule of mixtures (
ρ
eff
=
ρ
c
V
c
+
ρ
m
V
m
). The effective heat conductivity κ
eff and the coefficient of thermal expansion α
eff are given by
(21)
κ
eff
-
κ
m
κ
c
-
κ
m
=
V
c
1
+
V
m
(
(
κ
c
-
κ
m
)
/
3
κ
m
)
,
α
eff
-
α
m
α
c
-
α
m
=
(
1
/
K
eff
-
1
/
K
m
)
(
1
/
K
c
-
1
/
K
m
)
.Temperature-Dependent Material Property. The material properties that are temperature dependent are written as [41]
(22)
P
=
P
o
(
P
-
1
T
-
1
+
1
+
P
1
T
+
P
2
T
2
+
P
3
T
3
)
,
where P
o
,
P
-
1
,
P
1
,
P
2, and P
3 are the coefficients of temperature T and are unique to each constituent material phase.Temperature Distribution through the Thickness. The temperature variation is assumed to occur in the thickness direction only, and the temperature field is considered to be constant in the x
y-plane. In such a case, the temperature distribution along the thickness can be obtained by solving a steady state heat transfer problem:(23)
-
d
d
z
[
κ
(
z
)
d
T
d
z
]
=
0
,
T
=
{
T
c
at
z
=
h
2
,
T
m
at
z
=
-
h
2
.
The solution of (23) is obtained by means of a polynomial series [42] as
(24)
T
(
z
)
=
T
m
+
(
T
c
-
T
m
)
η
(
z
,
h
)
,
where
(25)
η
(
z
,
h
)
=
1
C
[
(
2
z
+
h
2
h
)
-
κ
c
m
(
n
+
1
)
κ
m
(
2
z
+
h
2
h
)
n
+
1
kkl
l
+
κ
c
m
2
(
2
n
+
1
)
κ
m
2
(
2
z
+
h
2
h
)
2
n
+
1
kkl
l
-
κ
c
m
3
(
3
n
+
1
)
κ
m
3
(
2
z
+
h
2
h
)
3
n
+
1
kkl
l
+
κ
c
m
4
(
4
n
+
1
)
κ
m
4
(
2
z
+
h
2
h
)
4
n
+
1
kkl
l
-
κ
c
m
5
(
5
n
+
1
)
κ
m
5
(
2
z
+
h
2
h
)
5
n
+
1
]
;
(26)
C
=
1
-
κ
c
m
(
n
+
1
)
κ
m
+
κ
c
m
2
(
2
n
+
1
)
κ
m
2
-
κ
c
m
3
(
3
n
+
1
)
κ
m
3
+
κ
c
m
4
(
4
n
+
1
)
κ
m
4
-
κ
c
m
5
(
5
n
+
1
)
κ
m
5
,
where κ
c
m
=
κ
c
-
κ
m.


## 2.1. Mindlin-Reissner Plate Theory

The Mindlin-Reissner plate theory, also known as the first-order shear deformation theory (FSDT), takes into account the shear deformation through the thickness. Using the Mindlin formulation, the displacementsu, v, and w at a point (
x
,
y
,
z
) in the plate (see Figure 1) from the medium surface are expressed as functions of the midplane displacements u
o, v
o, and w
o and independent rotations θ
x and θ
y of the normal in y
z and x
z planes, respectively, as follows:
(1)
u
(
x
,
y
,
z
,
t
)
=
u
o
(
x
,
y
,
t
)
+
z
θ
x
(
x
,
y
,
t
)
,
v
(
x
,
y
,
z
,
t
)
=
v
o
(
x
,
y
,
t
)
+
z
θ
y
(
x
,
y
,
t
)
,
w
(
x
,
y
,
z
,
t
)
=
w
o
(
x
,
y
,
t
)
,
where t is the time. The strains in terms of midplane deformation can be written as follows:
(2)
ε
=
{
ε
p
0
}
+
{
z
ε
b
ε
s
}
.
The midplane strains ε
p, the bending strain ε
b, and the shear strain ε
s in (2) are written as follows:
(3)
ε
p
=
{
u
o
,
x
v
o
,
y
u
o
,
y
+
v
o
,
x
}
,
ε
b
=
{
θ
x
,
x
θ
y
,
y
θ
x
,
y
+
θ
y
,
x
}
,
ε
s
=
{
θ
x
+
w
o
,
x
θ
y
+
w
o
,
y
}
,
where the subscript “comma” represents the partial derivative with respect to the spatial coordinate succeeding it. The membrane stress resultants N and the bending stress resultants M can be related to the membrane strains ε
p and bending strains ε
b through the following constitutive relations:
(4)
N
=
{
N
x
x
N
y
y
N
x
y
}
=
A
ε
p
+
B
ε
b
-
N
th
,
M
=
{
M
x
x
M
y
y
M
x
y
}
=
B
ε
p
+
D
b
ε
b
-
M
th
,
where the matrices A
=
A
i
j
,
B
=
B
i
j, and D
b
=
D
i
j
;
(
i
,
j
=
1,2
,
6
) are the extensional and bending-extensional coupling and bending stiffness coefficients and are defined as follows:
(5)
{
A
i
j
,
B
i
j
,
D
i
j
}
=
∫
-
h
/
2
h
/
2
‍
Q
¯
i
j
{
1
,
z
,
z
2
}
d
z
.
Similarly, the transverse shear force Q
=
{
Q
x
z
,
Q
y
z
} is related to the transverse shear strains ε
s through the following equation:
(6)
Q
i
j
=
E
i
j
ε
s
,
where E
=
E
i
j
=
∫
-
h
/
2
h
/
2
‍
Q
¯
i
j
υ
i
υ
j
d
z
;
(
i
,
j
=
4,5
) are the transverse shear stiffness coefficients, and υ
i, υ
j is the transverse shear coefficient for nonuniform shear strain distribution through the plate thickness. The stiffness coefficients Q
¯
i
j are defined as
(7)
Q
¯
11

=

Q
¯
22
=
E
(
z
)
1
-
ν
2
;
Q
¯
12
=
ν
E
(
z
)
1
-
ν
2
;
Q
¯
16
=
Q
¯
26
=
0
,
Q
¯
44
=
Q
¯
55
=
Q
¯
66
=
E
(
z
)
2
(
1
+
ν
)
,
where the modulus of elasticity E
(
z
) and Poisson’s ratio ν are given by (20). The thermal stress resultant N
th and the moment resultant M
th are
(8)
N
th
=
{
N
x
x
th
N
y
y
th
N
x
y
th
}
=
∫
-
h
/
2
h
/
2
Q
¯
i
j
α
(
z
,
T
)
{
1
1
0
}
△
T
(
z
)
d
z
,
M
th
=
{
M
x
x
th
M
y
y
th
M
x
y
th
}
=
∫
-
h
/
2
h
/
2
Q
¯
i
j
α
(
z
,
T
)
{
1
1
0
}
△
T
(
z
)
z
d
z
,
where the thermal coefficient of expansion α
(
z
,
T
) is given by (21), △
T
(
z
)
=
T
(
z
)
-
T
o is the temperature rise from the reference temperature, and T
o is the temperature at which there are no thermal strains. The strain energy function U is given by
(9)
U
(
δ
)
=
1
2
∫
Ω
‍
{
ε
p
T
A
ε
p
+
ε
p
T
B
ε
b
+
ε
b
T
B
ε
p
k
k
k
k
k
k
l
l
+
ε
b
T
D
ε
b
+
ε
s
T
E
ε
s
-
ε
b
T
N
th
-
ε
b
T
M
th
}
d
Ω
,
where δ
=
{
u
,
v
,
w
,
θ
x
,
θ
y
} is the vector of the degree of freedom associated with the displacement field in a finite element discretization. Following the procedure given in [40], the strain energy function U given in (9) can be rewritten as
(10)
U
(
δ
)
=
1
2
δ
T
K
δ
,
where K is the linear stiffness matrix. The kinetic energy of the plate is given by
(11)
T
(
δ
)
=
1
2
∫
Ω
‍
{
p
(
u
˙
o
2
+
v
˙
o
2
+
w
˙
o
2
)
+
I
(
θ
˙
x
2
+
θ
˙
y
2
)
}
d
Ω
,
where =
∫
-
h
/
2
h
/
2
‍
ρ
(
z
)
d
z
,
I
=
∫
-
h
/
2
h
/
2
‍
z
2
ρ
(
z
)
d
z, and ρ
(
z
) is the mass density that varies through the thickness of the plate. When the plate is subjected to a temperature field, this in turn results in in-plane stress resultants, N
th. The external work due to the in-plane stress resultants developed in the plate under a thermal load is given by
(12)
V
(
δ
)
=
∫
Ω
{
1
2
[
N
x
x
th
w
,
x
2
+
N
y
y
th
w
,
y
2
+
2
N
x
y
th
w
,
x
w
,
y
]
kkkl
+
h
2
24
[
N
x
x
th
(
θ
x
,
x
2
+
θ
y
,
x
2
)
+
N
y
y
2
(
θ
x
,
y
2
+
θ
y
,
y
2
)
kkkl
hh
+
2
N
x
y
th
(
θ
x
,
x
θ
x
,
y
+
θ
y
,
x
θ
y
,
y
)
]
1
2
}
d
Ω
.
Substituting (9)–(12) in Lagrange’s equation of motion, one obtains the following finite element equations:Static bending:(13)
K
δ
=
F
.
Free vibration:(14)
M
δ
¨
+
(
K
+
K
G
)
δ
=
0
.
Buckling analysis:
Mechanical buckling (prebuckling deformations are assumed to be zero or negligible in the analysis (including those coming from in-plane and out-of-plane coupling related to FGM and temperature variation through the thickness of the plate)):(15)
(
K
+
λ
M
K
G
)
δ
=
0
.
Thermal buckling:(16)
(
K
+
λ
T
K
G
)
δ
=
0
, where λ
M includes the critical value applied to in-plane mechanical loading, λ
T is the critical temperature difference, and K and K
G are the linear stiffness and geometric stiffness matrices, respectively. The critical temperature difference is computed using a standard eigenvalue algorithm.(a) Coordinate system of a rectangular FGM plate and (b) coordinate system of a skew plate.
(a)
(b)


## 2.2. Functionally Graded Material

A rectangular plate made of a mixture of ceramic and metal is considered with the coordinatesx
,
y along the in-plane directions and z along the thickness direction (see Figure 1). The material on the top surface (
z
=
h
/
2
) of the plate is ceramic rich and is graded to metal at the bottom surface of the plate (
z
=
-
h
/
2
) by a power law distribution. The effective properties of the FGM plate can be computed by using the rule of mixtures or by employing the Mori-Tanaka homogenization scheme. Let V
i
(
i
=
c
,
m
) be the volume fraction of the phase material. The subscripts c and m refer to ceramic and metal phases, respectively. The volume fraction of ceramic and metal phases is related by V
c
+
V
m
=
1, and V
c is expressed as
(17)
V
c
(
z
)
=
(
2
z
+
h
2
h
)
n
,
where n is the volume fraction exponent (
n
≥
0
), also known as the gradient index. The variation of the composition of ceramic and metal is linear for n
=
1, the value of n
=
0 represents a fully ceramic plate, and any other value of n yields a composite material with a smooth transition from ceramic to metal.Mori-Tanaka Homogenization Method. Based on the Mori-Tanaka homogenization method, the effective Young’s modulus and Poisson’s ratio are computed from the effective bulk modulus K and the effective shear modulus G as [41]
(18)
K
eff
-
K
m
K
c
-
K
m
=
V
c
1
+
V
m
(
3
(
K
c
-
K
m
)
/
(
3
K
m
+
4
G
m
)
)
,
G
eff
-
G
m
G
c
-
G
m
=
V
c
1
+
V
m
(
(
G
c
-
G
m
)
/
(
G
m
+
f
1
)
)
,
where
(19)
f
1
=
G
m
(
9
K
m
+
8
G
m
)
6
(
K
m
+
2
G
m
)
.
The effective Young’s modulus E
eff and Poisson’s ratio ν
eff can be computed from the following relations:
(20)
E
eff
=
9
K
eff
G
eff
3
K
eff
+
G
eff
,
ν
eff
=
3
K
eff
-
2
G
eff
2
(
3
K
eff
+
G
eff
)
.
The effective mass density ρ
eff is computed using the rule of mixtures (
ρ
eff
=
ρ
c
V
c
+
ρ
m
V
m
). The effective heat conductivity κ
eff and the coefficient of thermal expansion α
eff are given by
(21)
κ
eff
-
κ
m
κ
c
-
κ
m
=
V
c
1
+
V
m
(
(
κ
c
-
κ
m
)
/
3
κ
m
)
,
α
eff
-
α
m
α
c
-
α
m
=
(
1
/
K
eff
-
1
/
K
m
)
(
1
/
K
c
-
1
/
K
m
)
.Temperature-Dependent Material Property. The material properties that are temperature dependent are written as [41]
(22)
P
=
P
o
(
P
-
1
T
-
1
+
1
+
P
1
T
+
P
2
T
2
+
P
3
T
3
)
,
where P
o
,
P
-
1
,
P
1
,
P
2, and P
3 are the coefficients of temperature T and are unique to each constituent material phase.Temperature Distribution through the Thickness. The temperature variation is assumed to occur in the thickness direction only, and the temperature field is considered to be constant in the x
y-plane. In such a case, the temperature distribution along the thickness can be obtained by solving a steady state heat transfer problem:(23)
-
d
d
z
[
κ
(
z
)
d
T
d
z
]
=
0
,
T
=
{
T
c
at
z
=
h
2
,
T
m
at
z
=
-
h
2
.
The solution of (23) is obtained by means of a polynomial series [42] as
(24)
T
(
z
)
=
T
m
+
(
T
c
-
T
m
)
η
(
z
,
h
)
,
where
(25)
η
(
z
,
h
)
=
1
C
[
(
2
z
+
h
2
h
)
-
κ
c
m
(
n
+
1
)
κ
m
(
2
z
+
h
2
h
)
n
+
1
kkl
l
+
κ
c
m
2
(
2
n
+
1
)
κ
m
2
(
2
z
+
h
2
h
)
2
n
+
1
kkl
l
-
κ
c
m
3
(
3
n
+
1
)
κ
m
3
(
2
z
+
h
2
h
)
3
n
+
1
kkl
l
+
κ
c
m
4
(
4
n
+
1
)
κ
m
4
(
2
z
+
h
2
h
)
4
n
+
1
kkl
l
-
κ
c
m
5
(
5
n
+
1
)
κ
m
5
(
2
z
+
h
2
h
)
5
n
+
1
]
;
(26)
C
=
1
-
κ
c
m
(
n
+
1
)
κ
m
+
κ
c
m
2
(
2
n
+
1
)
κ
m
2
-
κ
c
m
3
(
3
n
+
1
)
κ
m
3
+
κ
c
m
4
(
4
n
+
1
)
κ
m
4
-
κ
c
m
5
(
5
n
+
1
)
κ
m
5
,
where κ
c
m
=
κ
c
-
κ
m.


## 3. Cell-Based Smoothed Finite Element Method with Discrete Shear Gap Technique

In this study, three-noded triangular element with five degrees of freedom (DOFs)δ
=
{
u
,
v
,
w
,
θ
x
,
θ
y
} is employed. The displacement is approximated by
(27)
u
h
=
∑
I
‍
N
I
δ
I
,
where δ
I are the nodal DOFs and N
I are the standard finite element shape functions given by
(28)
N
=
[
1
-
ξ
-
η
,
η
,
ξ
]
.In this work, the cell-based smoothed finite element method (CSFEM) is combined with stabilized discrete shear gap method (DSG) for three-noded triangular element, called “cell-based discrete shear gap method (CS-DSG3)” [39]. The cell-based smoothing technique decreases the computational complexity, whilst DSG suppresses the shear locking phenomenon when the present formulation is applied to thin plates. Interested readers are referred to the literature and references therein for the description of cell-based smoothing technique [24, 26] and DSG method [43]. In the CS-DSG3, each triangular element is divided into three subtriangles. The displacement vector at the center node is assumed to be the simple average of the three displacement vectors of the three field nodes. In each subtriangle, the stabilized DSG3 is used to compute the strains and also to avoid the transverse shear locking. Then the strain smoothing technique on the whole triangular element is used to smooth the strains on the three subtriangles. Consider a typical triangular element Ω
e as shown in Figure 2. This is first divided into three subtriangles △
1
,
△
2, and △
3 such that Ω
e
=
⋃
i
=
1
3
△
i. The coordinates of the center point x
o
=
(
x
o
,
y
o
) are given by
(29)
(
x
o
,
y
o
)
=
1
3
(
x
I
,
y
I
)
.
The displacement vector of the center point is assumed to be a simple average of the nodal displacements as
(30)
δ
e
O
=
1
3
δ
e
I
.
The constant membrane strains, the bending strains, and the shear strains for subtriangle △
1 are given by
(31)
ε
p
=
[
p
1
△
1
p
2
△
1
p
3
△
1
]
{
δ
e
O
δ
e
1
δ
e
2
}
,
ε
b
=
[
b
1
△
1
b
2
△
1
b
3
△
1
]
{
δ
e
O
δ
e
1
δ
e
2
}
,
ε
s
=
[
s
1
△
1
s
2
△
1
s
3
△
1
]
{
δ
e
O
δ
e
1
δ
e
2
}
.
Upon substituting the expression for δ
e
O in (31), we obtain
(32)
ε
p
△
1
=
[
1
3
p
1
△
1
+
p
2
△
1
1
3
p
1
△
1
+
p
3
△
1
1
3
p
1
△
1
]
×
{
δ
e
1
δ
e
2
δ
e
3
}
=
B
p
△
1
δ
e
,
ε
b
△
1
=
[
1
3
b
1
△
1
+
b
2
△
1
1
3
b
1
△
1
+
b
3
△
1
1
3
b
1
△
1
]
×
{
δ
e
1
δ
e
2
δ
e
3
}
=
B
b
△
1
δ
e
,
ε
s
△
1
=
[
1
3
s
1
△
1
+
s
2
△
1
1
3
s
1
△
1
+
s
3
△
1
1
3
s
1
△
1
]
×
{
δ
e
1
δ
e
2
δ
e
3
}
=
B
s
△
1
δ
e
,
where p
i
,
(
i
=
1
,
2
,
3
), b
i
,
(
i
=
1
,
2
,
3
), and s
i
,
(
i
=
1
,
2
,
3
) are given by(33)
B
p
=
1
2
A
e
[
b
-
c
0
0
0
0
0
d
-
a
0
0
0
d
-
a
b
-
c
0
0
0
︸
p
1
c
0
0
0
0
0
-
d
0
0
0
-
d
c
0
0
0
︸
p
2
-
b
0
0
0
0
a
0
0
0
0
a
-
b
0
0
0
︸
p
3
]
,
B
b
=
1
2
A
e
[
0
0
0
b
-
c
0
0
0
0
0
d
-
a
0
0
0
d
-
a
b
-
c
︸
b
1
0
0
0
c
0
0
0
0
0
-
d
0
0
0
-
d
c
︸
b
2
0
0
0
-
b
0
0
0
0
0
a
0
0
0
a
-
b
︸
b
3
]
,
B
s
=
1
2
A
e
[
0
0
b
-
c
A
e
0
0
0
d
-
a
0
A
e
︸
s
1
0
0
c
a
c
/
2
b
c
/
2
0
0
-
d
-
a
d
/
2
-
b
d
/
2
︸
s
2
0
0
-
b
-
b
d
/
2
-
b
c
/
2
0
0
a
a
d
/
2
a
c
/
2
︸
s
3
]
,where a
=
x
2
-
x
1
;
b
=
y
2
-
y
1
;
c
=
y
3
-
y
1
;
d
=
x
3
-
x
1 (see Figure 3). A
e is the area of the triangular element and B
s is the altered shear strains [43]. The strain-displacement matrix for the other two triangles can be obtained by cyclic permutation. Now applying the cell-based strain smoothing [26], the constant membrane strains, the bending strains, and the shear strains are, respectively, employed to create a smoothed membrane strain ε
¯
p, smoothed bending strain ε
¯
b, and smoothed shear strain ε
¯
s on the triangular element Ω
e as follows:
(34)
ε
¯
p
=
∫
Ω
e
ε
b
Φ
e
(
x
)
d
Ω
=
∑
i
=
1
3
‍
ε
p
△
i
∫
△
i
Φ
e
(
x
)
d
Ω
,
ε
¯
b
=
∫
Ω
e
ε
b
Φ
e
(
x
)
d
Ω
=
∑
i
=
1
3
‍
ε
b
△
i
∫
△
i
Φ
e
(
x
)
d
Ω
,
ε
¯
s
=
∫
Ω
e
ε
s
Φ
e
(
x
)
d
Ω
=
∑
i
=
1
3
‍
ε
s
△
i
∫
△
i
Φ
e
(
x
)
d
Ω
,
where Φ
e
(
x
) is a given smoothing function that satisfies
(35)
Φ
(
x
)
=
{
1
A
c
x
∈
Ω
c
0
x
∉
Ω
c
,
where A
c is the area of the triangular element. The smoothed membrane strain, the smoothed bending strain, and the smoothed shear strain are then given by
(36)
{
ε
¯
p
,
ε
¯
b
,
ε
¯
s
}
=
∑
i
=
1
3
A
△
i
{
ε
p
△
i
,
ε
b
△
i
,
ε
s
△
i
}
A
e
.
The smoothed elemental stiffness matrix is given by
(37)
K
=
∫
Ω
e
B
¯
p
A
B
¯
p
T
+
B
¯
p
B
B
¯
b
T
+
B
¯
b
B
B
¯
p
T
+
B
¯
b
D
B
¯
b
T
+
B
¯
s
E
B
¯
s
T
d
Ω
=
(
B
¯
p
A
B
¯
p
T
+
B
¯
p
B
B
¯
b
T
+
B
¯
b
B
B
¯
p
T
+
B
¯
b
D
B
¯
b
T
+
B
¯
s
E
B
¯
s
T
)
A
e
,
where B
¯
p
,
B
¯
b, and B
¯
s are the smoothed strain-displacement matrix. The mass matrix M is computed by following the conventional finite element procedure. To further improve the accuracy of the solution and to stabilize the shear force oscillation, the shear stiffness coefficients are multiplied by the following factor:
(38)
ShearFac
=
h
3
h
2
+
α
h
e
2
,
where α is a positive constant and h
e is the longest length of the edge of an element.Figure 2
A triangular element is divided into three subtriangles.△
1
,
△
2, and △
3 are the subtriangles created by connecting the central point O with three field nodes.Figure 3
Three-noded triangular element and local coordinates in discrete shear gap method.


## 4. Numerical Examples

In this section, we present the static bending response, the linear free vibration, and the buckling analysis of FGM plates using cell-based smoothed finite element method with discrete shear gap technique. The effect of various parameters, namely, material gradient indexn, skewness of the plate ψ, the plate aspect ratio a
/
b, the plate thickness a
/
h, and boundary conditions on the global response, is numerically studied. The top surface of the plate is ceramic rich and the bottom surface of the plate is metal rich. Here, the modified shear correction factor obtained based on energy equivalence principle as outlined in [44] is used. The boundary conditions for simply supported and clamped cases are as follows:simply supported boundary condition:(39)
u
o
=
w
o
=
θ
y
=
0
on
⁡
x
=
0
,
a
;
v
o
=
w
o
=
θ
x
=
0
on
⁡
y
=
0
,
b
;
clamped boundary condition:(40)
u
o
=
w
o
=
θ
y
=
v
o
=
θ
x
=
0
on
⁡
x
=
0
,
a
,
y
=
0
,
b
.Skew Boundary Transformation. For skew plates, the edges of the boundary elements may not be parallel to the global axes (
x
,
y
,
z
). In order to specify the boundary conditions on skew edges, it is necessary to use the edge displacements (
u
o
′
,
v
o
′
,
w
o
′
) and so forth in a local coordinate system (
x
′
,
y
′
,
z
′
) (see Figure 1). The element matrices corresponding to the skew edges are transformed from global axes to local axes on which the boundary conditions can be conveniently specified. The relation between the global and the local degrees of freedom of a particular node is obtained by
(41)
δ
=
L
g
δ
′
,
where δ and δ
′ are the generalized displacement vector in the global and the local coordinate system, respectively. The nodal transformation matrix for a node I on the skew boundary is given by
(42)
L
g
=
[
cos
⁡
ψ
sin
ψ
0
0
0
-
sin
ψ
cos
⁡
ψ
0
0
0
0
0
1
0
0
0
0
0
cos
⁡
ψ
sin
ψ
0
0
0
-
sin
ψ
cos
⁡
ψ
]
,
where ψ defines the skewness of the plate.


### 4.1. Static Bending

Let us consider aAl
/
ZrO
2 FGM square plate with length-to-thickness a
/
h
= 5, subjected to a uniform load with fully simply supported (SSSS) boundary conditions. The Young’s modulus for ZrO
2 is E
c
= 151 GPa and for aluminum is E
m
= 70 GPa. Poisson’s ratio is chosen as constant, ν
= 0.3. Table 1 compares the results from the present formulation with other approaches available in the literature [22, 45–47] and a very good agreement can be observed. Next, we illustrate the performance of the present formulation for thin plate problems. A simply supported square plate subjected to uniform load is considered, while the length-to-thickness (
a
/
h
) varies from 5 to 10
4. Three individual approaches are employed: discrete shear gap method referred to as DSG3, the cell-based smoothed finite element method with discrete shear gap technique (CSDSG3), and the stabilized CSDSG3 (where the shear stiffness coefficients are multiplied by the stabilization factor). The normalized center deflection w
¯
c
=
100
w
c
(
E
m
h
3
/
12
(
1
-
ν
2
)
p
a
4
) is shown in Figure 4. It is observed that the DSG3 results are subjected to shear locking when the plate becomes thin (
a
/
h
>
100
). However, the present formulation, CSDSG3 with stabilization, is less sensitive to shear locking.Table 1
The normalized center deflectionw
¯
c
=
100
w
c
(
E
c
h
3
/
12
(
1
-
ν
2
)
p
a
4
) for a simply supported Al/ZrO2-1 FGM square plate with a
/
h
= 5, subjected to a uniformly distributed load p.
Method
Gradient index,n
0
1
2
4
×
4
0.1443
0.2356
0.2644
8
×
8
0.1648
0.2703
0.3029
16
×
16
0.1701
0.2795
0.3131
32
×
32
0.1714
0.2819
0.3158
40
×
40
0.1716
0.2822
0.3161
NS-DSG3 [47]
0.1721
0.2716
0.3107
ES-DSG3 [47]
0.1700
0.2680
0.3066
MLPG [45]
0.1671
0.2905
0.3280
k
p-Ritz [46]
0.1722
0.2811
0.3221
MITC4 [47]
0.1715
0.2704
0.3093
IGA-Quadratic [22]
0.1717
0.2719
0.3115Figure 4
The normalized center deflection as a function of normalized plate thickness for a simply supported square FGM plate subjected to a uniform load.


### 4.2. Free Flexural Vibrations

In this section, the free flexural vibration characteristics of FGM plates with and without centrally located cutout in thermal environment are studied numerically. Figure5 shows the geometry of the plate with a centrally located circular cutout. In all cases, we present the nondimensionalized free flexural frequency defined as follows, unless otherwise stated:
(43)
ω
¯
=
ω
a
2
ρ
c
h
D
c
,
where ω is the natural frequency and ρ
c
,
D
c
=
E
c
h
3
/
12
(
1
-
ν
2
) are the mass density and the flexural rigidity of the ceramic phase. The FGM plate considered here is made up of silicon nitride (Si
3
N
4) and stainless steel (SUS304). The material is considered to be temperature dependent, and the temperature coefficients corresponding to Si
3
N
4/SUS304 are listed in Table 2 [13, 41]. The mass density (
ρ
) and the thermal conductivity (
κ
) are ρ
c
= 2370 kg/m
3and κ
c
= 9.19 W/mK for Si
3
N
4 and ρ
m
= 8166 kg/m
3and κ
m
= 12.04 W/mK for SUS304. Poisson’s ratio ν is assumed to be constant and taken as 0.28 for the current study [41]. Before proceeding with a detailed study on the effect of gradient index on the natural frequencies, the formulation developed herein is validated against available analytical/numerical solutions pertaining to the linear frequencies of a FGM plate in thermal environment and a FGM plate with a centrally located circular cutout. The computed frequencies (a) for a square simply supported FGM plate in thermal environment with a
/
h
= 10 are given in Table 3 and (b) the mesh convergence and comparison of linear frequencies for a square plate with circular cutout are given in Tables 4 and 5. It can be seen that the numerical results from the present formulation are found to be in very good agreement with the existing solutions. For the uniform temperature case, the material properties are evaluated at T
c
=
T
m
= 300 K. The temperature is assumed to vary only in the thickness direction and is determined by (24). The temperature for the ceramic surface is varied, whilst a constant value on the metallic surface is maintained (
T
m
= 300 K) to subject a thermal gradient. The geometric stiffness matrix is computed from the in-plane stress resultants due to the applied thermal gradient. The geometric stiffness matrix is then added to the stiffness matrix, and the eigenvalue problem is solved. The effect of the material gradient index is also shown in Tables 3 and 5 and the influence of a centrally located cutout is shown in Tables 4 and 5. The combined effect of increasing the temperature and the gradient index is to lower the fundamental frequency and this is due to the increase in the metallic volume fraction. Figure 6 shows the influence of the cutout size on the frequency for a plate in thermal environment (
△
T
= 100 K). The frequency increases with increasing cutout size. This can be attributed to the decrease in stiffness due to the presence of the cutout. Also, it can be seen that with increasing gradient index, the frequency decreases. In this case, the decrease in the frequency is due to the increase in the metallic volume fraction. It is observed that the combined effect of increasing the gradient index and the cutout size is to lower the fundamental frequency. Increasing the thermal gradient further decreases the fundamental frequency.Table 2
Temperature-dependent coefficient for material Si3N4/SUS304 [13, 41].
Material
Property
P
o
P
-
1
P
1
P
2
P
3
Si3N4
E (Pa)
348.43
e
9
0.0
-
3.070
e
-
4
2.160
e
-
7
-
8.946
e
-
11
α (1/K)
5.8723
e
-
6
0.0
9.095
e
-
4
0.0
0.0
SUS304
E (Pa)
201.04
e
9
0.0
3.079
e
-
4
-
6.534
e
-
7
0.0
α (1/K)
12.330
e
-
6
0.0
8.086
e
-
4
0.0
0.0Table 3
The first normalized frequency parameterω
¯ for a fully simply supported Si3N4/SUS304 FGM square plate with a
/
h
= 10 in thermal environment.
T
c, T
m
Gradient index,n
0
1
5
10
300, 300
Present
18.3570
11.0690
9.0260
8.5880
[48]
18.3731
11.0288
9.0128
8.5870
400, 300
Present
17.9778
10.7979
8.8626
8.3182
[48]
17.9620
10.7860
8.7530
8.3090
600, 300
Present
17.1205
10.1679
8.1253
7.6516
[48]
17.1050
10.1550
8.1150
7.6420Table 4
Convergence of fundamental frequency(
Ω
=
[
ω
2
ρ
c
h
a
4
/
D
c
(
1
-
ν
2
)
]
1
/
4
) with mesh size for an isotropic plate with a central cutout.
Number of nodes
Mode 1
Mode 2
333
6.1025
8.6297
480
6.0805
8.5595
719
6.0663
8.5192
1271
6.0560
8.4852
[49]
6.1725
8.6443
[50]
6.2110
8.7310Table 5
Comparison of fundamental frequency for a simply supported FGM plate witha
/
h
= 5 and r
/
a
= 0.2.
T
c
Gradient index,n
0
1
2
5
10
300
[49]
17.6855
10.6681
9.6040
8.7113
8.2850
Present
17.7122
10.6845
9.6188
8.7246
8.2976
400
[37]
17.4690
10.5174
9.4618
8.5738
8.1484
Present
17.5488
10.5775
9.5197
8.6309
8.2059Figure 5
Plate with a centrally located circular cutout.r is the radius of the circular cutout.Figure 6
Effect of the cutout size on the fundamental frequency(
Ω
) for a square simply supported FGM plate with a central circular cutout in thermal environment △
T
= 100 K (
T
c
= 400 K, T
m
= 300 K) for different gradient index n.


### 4.3. Buckling Analysis

In this section, we present the mechanical and thermal buckling behaviour of functionally graded skew plates.Mechanical Buckling. The FGM plate considered here consists of aluminum (Al) and zirconium dioxide (ZrO2). The material is considered to be temperature independent. The Young’s modulus (
E
) forZr
O
2 is E
c
= 151 GPa and for Al is E
m
= 70 GPa. For mechanical buckling, we consider both uni- and biaxial mechanical loads on the FGM plates. In all cases, we present the critical buckling parameters as follows, unless otherwise specified:
(44)
λ
cru
=
N
x
x
cr
0
b
2
π
2
D
c
,
λ
crb
=
N
y
y
cr
0
b
2
π
2
D
c
,
where λ
cru and λ
crb are the critical buckling parameters for uni- and biaxial load, respectively; D
c
=
E
c
h
3
/
(
12
(
1
-
ν
2
)
). The critical buckling loads evaluated by varying the skew angle of the plate and volume fraction index and considering mechanical loads (uni- and biaxial compressive loads) are shown in Table 6 for a
/
h
= 100. The efficacy of the present formulation is demonstrated by comparing our results with those in [14]. It can be seen that increasing the gradient index decreases the critical buckling load. A very good agreement in the results can be observed. It is also observed that the decrease in the critical value is significant for the material gradient index n
≤
2 and that further increase in n yields less reduction in the critical value, irrespective of the skew angle. The effect of the plate aspect ratio and the gradient index on the critical buckling load is shown in Figure 7 for a simply supported FGM plate under uniaxial mechanical load. It is observed that the combined effect of increasing the gradient index and the plate aspect ratio is to lower the critical buckling load. Table 7 presents the critical buckling parameter for a simply supported FGM with a centrally located circular cutout with r
/
a
= 0.2. It can be seen that the present formulation yields comparable results. The effect of increasing the gradient index is to lower the critical buckling load. This can be attributed to the stiffness degradation due to increase in the metallic volume fraction. Figure 8 shows the influence of a centrally located circular cutout and the gradient index on the critical buckling load under two different boundary conditions, namely, all edges simply supported and all edges clamped. In this case, the plate is subjected to a uniaxial compressive load. It can be seen that increasing the gradient index decreases the critical buckling load due to increasing metallic volume fraction, whilst increasing the cutout radius decreases the critical buckling load in the case of simply supported boundary conditions. This can be attributed to the stiffness degradation due to the presence of a cutout and the simply supported boundary condition. In case of the clamped boundary condition, the critical buckling load first decreases with increasing cutout radius due to stiffness degradation. Upon further increase, the critical buckling load increases. This is because the clamped boundary condition adds stiffness to the system which overcomes the stiffness reduction due to the presence of a cutout.Table 6
Critical buckling parameters for a thin simply supported FGM skew plate witha
/
h
= 100 and a
/
b
= 1.
Skew angle
λ
cr
Gradient index,n
0
1
5
10
[14]
Present
[43]
Present
0°
λ
cru
4.0010
4.0034
1.7956
1.8052
1.2624
1.0846
λ
crb
2.0002
2.0017
0.8980
0.9028
0.6312
0.5423
15°
λ
cru
4.3946
4.4007
1.9716
1.9799
1.3859
1.1915
λ
crb
2.1154
2.1187
0.9517
0.9561
0.6683
0.5741
30°
λ
cru
5.8966
5.9317
2.6496
2.6496
1.8586
1.6020
λ
crb
2.5365
2.5491
1.1519
1.1520
0.8047
0.6909Table 7
Comparison of critical buckling loadλ
cru
=
N
x
x
cr
o
b
2
/
π
2
D
m for a simply supported FGM plate with a
/
h
= 100 and r
/
a
= 0.2. The effective material properties are computed by rule of mixtures. In order to be consistent with the literature, the properties of the metallic phase are used for normalization.
Gradient index,n
[51]
Present
% difference
0
5.2611
5.2831
−0.42
0.2
4.6564
4.6919
−0.76
1
3.6761
3.663
0.36
2
3.3672
3.3961
−0.86
5
3.1238
3.1073
0.53
10
2.9366
2.8947
1.43Figure 7
Effect of plate aspect ratioa
/
b and gradient index on the critical buckling load for a simply supported FGM plate under uniaxial compression with a
/
h
=
10.Figure 8
Variation of the critical buckling load,λ
cru
=
N
x
x
cr
o
b
2
/
π
2
D
c, with cutout dimensions for a square FGM plate with a
/
h
=
10 subjected to uniaxial compressive loading for different gradient index n and various boundary conditions.Thermal Buckling. The thermal buckling behaviour of simply supported functionally graded skew plate is studied next. The top surface is ceramic rich and the bottom surface is metal rich. The FGM plate considered here consists of aluminum and alumina. Young’s modulus, the thermal conductivity, and the coefficient of thermal expansion for alumina are E
c
= 380 GPa, κ
c
=10.4 W/mK, and α
c
=
7.4
×
10
-
6 1/°C and for aluminum E
m
= 70 GPa, κ
m
= 204 W/mK, and α
m
=
23
×
10
-
6 1/°C, respectively. Poisson’s ratio is chosen as constant, ν
= 0.3. The temperature rise of T
m
= 5°C in the metal-rich surface of the plate is assumed in the present study. In addition to nonlinear temperature distribution across the plate thickness, the linear case is also considered in the present analysis by truncating the higher order terms in (25). The plate is of uniform thickness and simply supported on all four edges. Table 8 shows the convergence of the critical buckling temperature with mesh size for different gradient index, n. It can be seen that the results from the present formulation are in very good agreement with the available solution. The influence of the plate aspect ratio a
/
b and the skew angle ψ on the critical buckling temperature for a simply supported square FGM plate are shown in Figures 9 and 10. It is seen that increasing the plate aspect ratio decreases the critical buckling temperature for both linear and nonlinear temperature distribution through the thickness. The critical buckling temperature increases with increase in the skew angle. The influence of the gradient index n is also shown in Figure 10. It is seen that with increasing gradient index, n, the critical buckling temperature decreases. This is due to the increase in the metallic volume fraction that degrades the overall stiffness of the structure. Figure 11 shows the influence of the cutout radius and the material gradient index on the critical buckling temperature. Both linear and nonlinear temperature distribution through the thickness are assumed. Again, it is seen that the combined effect of increasing the gradient index n and the cutout radius r
/
a is to lower the buckling temperature. For gradient index n
= 0, there is no difference between the linear and the nonlinear temperature distribution through the thickness as the material is homogeneous through the thickness, while, for n
>
0, the material is heterogeneous through the thickness with different thermal property.Table 8
Convergence of the critical buckling temperature for a simply supported FGM skew plate witha
/
h
= 10 and a
/
b
= 1. Nonlinear temperature rise through the thickness of the plate is assumed.
Mesh
Gradient index,n
0
1
5
10
8
×
8
3383.40
2054.61
1539.24
1496.36
16
×
16
3286.90
1995.07
1495.25
1453.99
32
×
32
3263.91
1980.96
1484.76
1443.86
40
×
40
3261.17
1979.30
1483.51
1442.60
[9]
3257.47
1977.01
1481.83
1441.02Figure 9
Critical buckling temperature as a function of plate aspect ratioa
/
b with linear and nonlinear temperature distribution through the thickness with gradient index n
= 5.Figure 10
Critical buckling temperature as a function of skew angleψ for a simply supported square FGM plate with a
/
h
=
10. Both linear and nonlinear temperature distribution through the thickness are assumed.Figure 11
Influence of cutout size on the critical buckling temperature for a square simply supported FGM plate with a centrally located circular cutout witha
/
h
= 10 for various gradient indices n. Linear and nonlinear temperature distribution through the thickness are assumed.


## 4.1. Static Bending

Let us consider aAl
/
ZrO
2 FGM square plate with length-to-thickness a
/
h
= 5, subjected to a uniform load with fully simply supported (SSSS) boundary conditions. The Young’s modulus for ZrO
2 is E
c
= 151 GPa and for aluminum is E
m
= 70 GPa. Poisson’s ratio is chosen as constant, ν
= 0.3. Table 1 compares the results from the present formulation with other approaches available in the literature [22, 45–47] and a very good agreement can be observed. Next, we illustrate the performance of the present formulation for thin plate problems. A simply supported square plate subjected to uniform load is considered, while the length-to-thickness (
a
/
h
) varies from 5 to 10
4. Three individual approaches are employed: discrete shear gap method referred to as DSG3, the cell-based smoothed finite element method with discrete shear gap technique (CSDSG3), and the stabilized CSDSG3 (where the shear stiffness coefficients are multiplied by the stabilization factor). The normalized center deflection w
¯
c
=
100
w
c
(
E
m
h
3
/
12
(
1
-
ν
2
)
p
a
4
) is shown in Figure 4. It is observed that the DSG3 results are subjected to shear locking when the plate becomes thin (
a
/
h
>
100
). However, the present formulation, CSDSG3 with stabilization, is less sensitive to shear locking.Table 1
The normalized center deflectionw
¯
c
=
100
w
c
(
E
c
h
3
/
12
(
1
-
ν
2
)
p
a
4
) for a simply supported Al/ZrO2-1 FGM square plate with a
/
h
= 5, subjected to a uniformly distributed load p.
Method
Gradient index,n
0
1
2
4
×
4
0.1443
0.2356
0.2644
8
×
8
0.1648
0.2703
0.3029
16
×
16
0.1701
0.2795
0.3131
32
×
32
0.1714
0.2819
0.3158
40
×
40
0.1716
0.2822
0.3161
NS-DSG3 [47]
0.1721
0.2716
0.3107
ES-DSG3 [47]
0.1700
0.2680
0.3066
MLPG [45]
0.1671
0.2905
0.3280
k
p-Ritz [46]
0.1722
0.2811
0.3221
MITC4 [47]
0.1715
0.2704
0.3093
IGA-Quadratic [22]
0.1717
0.2719
0.3115Figure 4
The normalized center deflection as a function of normalized plate thickness for a simply supported square FGM plate subjected to a uniform load.


## 4.2. Free Flexural Vibrations

In this section, the free flexural vibration characteristics of FGM plates with and without centrally located cutout in thermal environment are studied numerically. Figure5 shows the geometry of the plate with a centrally located circular cutout. In all cases, we present the nondimensionalized free flexural frequency defined as follows, unless otherwise stated:
(43)
ω
¯
=
ω
a
2
ρ
c
h
D
c
,
where ω is the natural frequency and ρ
c
,
D
c
=
E
c
h
3
/
12
(
1
-
ν
2
) are the mass density and the flexural rigidity of the ceramic phase. The FGM plate considered here is made up of silicon nitride (Si
3
N
4) and stainless steel (SUS304). The material is considered to be temperature dependent, and the temperature coefficients corresponding to Si
3
N
4/SUS304 are listed in Table 2 [13, 41]. The mass density (
ρ
) and the thermal conductivity (
κ
) are ρ
c
= 2370 kg/m
3and κ
c
= 9.19 W/mK for Si
3
N
4 and ρ
m
= 8166 kg/m
3and κ
m
= 12.04 W/mK for SUS304. Poisson’s ratio ν is assumed to be constant and taken as 0.28 for the current study [41]. Before proceeding with a detailed study on the effect of gradient index on the natural frequencies, the formulation developed herein is validated against available analytical/numerical solutions pertaining to the linear frequencies of a FGM plate in thermal environment and a FGM plate with a centrally located circular cutout. The computed frequencies (a) for a square simply supported FGM plate in thermal environment with a
/
h
= 10 are given in Table 3 and (b) the mesh convergence and comparison of linear frequencies for a square plate with circular cutout are given in Tables 4 and 5. It can be seen that the numerical results from the present formulation are found to be in very good agreement with the existing solutions. For the uniform temperature case, the material properties are evaluated at T
c
=
T
m
= 300 K. The temperature is assumed to vary only in the thickness direction and is determined by (24). The temperature for the ceramic surface is varied, whilst a constant value on the metallic surface is maintained (
T
m
= 300 K) to subject a thermal gradient. The geometric stiffness matrix is computed from the in-plane stress resultants due to the applied thermal gradient. The geometric stiffness matrix is then added to the stiffness matrix, and the eigenvalue problem is solved. The effect of the material gradient index is also shown in Tables 3 and 5 and the influence of a centrally located cutout is shown in Tables 4 and 5. The combined effect of increasing the temperature and the gradient index is to lower the fundamental frequency and this is due to the increase in the metallic volume fraction. Figure 6 shows the influence of the cutout size on the frequency for a plate in thermal environment (
△
T
= 100 K). The frequency increases with increasing cutout size. This can be attributed to the decrease in stiffness due to the presence of the cutout. Also, it can be seen that with increasing gradient index, the frequency decreases. In this case, the decrease in the frequency is due to the increase in the metallic volume fraction. It is observed that the combined effect of increasing the gradient index and the cutout size is to lower the fundamental frequency. Increasing the thermal gradient further decreases the fundamental frequency.Table 2
Temperature-dependent coefficient for material Si3N4/SUS304 [13, 41].
Material
Property
P
o
P
-
1
P
1
P
2
P
3
Si3N4
E (Pa)
348.43
e
9
0.0
-
3.070
e
-
4
2.160
e
-
7
-
8.946
e
-
11
α (1/K)
5.8723
e
-
6
0.0
9.095
e
-
4
0.0
0.0
SUS304
E (Pa)
201.04
e
9
0.0
3.079
e
-
4
-
6.534
e
-
7
0.0
α (1/K)
12.330
e
-
6
0.0
8.086
e
-
4
0.0
0.0Table 3
The first normalized frequency parameterω
¯ for a fully simply supported Si3N4/SUS304 FGM square plate with a
/
h
= 10 in thermal environment.
T
c, T
m
Gradient index,n
0
1
5
10
300, 300
Present
18.3570
11.0690
9.0260
8.5880
[48]
18.3731
11.0288
9.0128
8.5870
400, 300
Present
17.9778
10.7979
8.8626
8.3182
[48]
17.9620
10.7860
8.7530
8.3090
600, 300
Present
17.1205
10.1679
8.1253
7.6516
[48]
17.1050
10.1550
8.1150
7.6420Table 4
Convergence of fundamental frequency(
Ω
=
[
ω
2
ρ
c
h
a
4
/
D
c
(
1
-
ν
2
)
]
1
/
4
) with mesh size for an isotropic plate with a central cutout.
Number of nodes
Mode 1
Mode 2
333
6.1025
8.6297
480
6.0805
8.5595
719
6.0663
8.5192
1271
6.0560
8.4852
[49]
6.1725
8.6443
[50]
6.2110
8.7310Table 5
Comparison of fundamental frequency for a simply supported FGM plate witha
/
h
= 5 and r
/
a
= 0.2.
T
c
Gradient index,n
0
1
2
5
10
300
[49]
17.6855
10.6681
9.6040
8.7113
8.2850
Present
17.7122
10.6845
9.6188
8.7246
8.2976
400
[37]
17.4690
10.5174
9.4618
8.5738
8.1484
Present
17.5488
10.5775
9.5197
8.6309
8.2059Figure 5
Plate with a centrally located circular cutout.r is the radius of the circular cutout.Figure 6
Effect of the cutout size on the fundamental frequency(
Ω
) for a square simply supported FGM plate with a central circular cutout in thermal environment △
T
= 100 K (
T
c
= 400 K, T
m
= 300 K) for different gradient index n.


## 4.3. Buckling Analysis

In this section, we present the mechanical and thermal buckling behaviour of functionally graded skew plates.Mechanical Buckling. The FGM plate considered here consists of aluminum (Al) and zirconium dioxide (ZrO2). The material is considered to be temperature independent. The Young’s modulus (
E
) forZr
O
2 is E
c
= 151 GPa and for Al is E
m
= 70 GPa. For mechanical buckling, we consider both uni- and biaxial mechanical loads on the FGM plates. In all cases, we present the critical buckling parameters as follows, unless otherwise specified:
(44)
λ
cru
=
N
x
x
cr
0
b
2
π
2
D
c
,
λ
crb
=
N
y
y
cr
0
b
2
π
2
D
c
,
where λ
cru and λ
crb are the critical buckling parameters for uni- and biaxial load, respectively; D
c
=
E
c
h
3
/
(
12
(
1
-
ν
2
)
). The critical buckling loads evaluated by varying the skew angle of the plate and volume fraction index and considering mechanical loads (uni- and biaxial compressive loads) are shown in Table 6 for a
/
h
= 100. The efficacy of the present formulation is demonstrated by comparing our results with those in [14]. It can be seen that increasing the gradient index decreases the critical buckling load. A very good agreement in the results can be observed. It is also observed that the decrease in the critical value is significant for the material gradient index n
≤
2 and that further increase in n yields less reduction in the critical value, irrespective of the skew angle. The effect of the plate aspect ratio and the gradient index on the critical buckling load is shown in Figure 7 for a simply supported FGM plate under uniaxial mechanical load. It is observed that the combined effect of increasing the gradient index and the plate aspect ratio is to lower the critical buckling load. Table 7 presents the critical buckling parameter for a simply supported FGM with a centrally located circular cutout with r
/
a
= 0.2. It can be seen that the present formulation yields comparable results. The effect of increasing the gradient index is to lower the critical buckling load. This can be attributed to the stiffness degradation due to increase in the metallic volume fraction. Figure 8 shows the influence of a centrally located circular cutout and the gradient index on the critical buckling load under two different boundary conditions, namely, all edges simply supported and all edges clamped. In this case, the plate is subjected to a uniaxial compressive load. It can be seen that increasing the gradient index decreases the critical buckling load due to increasing metallic volume fraction, whilst increasing the cutout radius decreases the critical buckling load in the case of simply supported boundary conditions. This can be attributed to the stiffness degradation due to the presence of a cutout and the simply supported boundary condition. In case of the clamped boundary condition, the critical buckling load first decreases with increasing cutout radius due to stiffness degradation. Upon further increase, the critical buckling load increases. This is because the clamped boundary condition adds stiffness to the system which overcomes the stiffness reduction due to the presence of a cutout.Table 6
Critical buckling parameters for a thin simply supported FGM skew plate witha
/
h
= 100 and a
/
b
= 1.
Skew angle
λ
cr
Gradient index,n
0
1
5
10
[14]
Present
[43]
Present
0°
λ
cru
4.0010
4.0034
1.7956
1.8052
1.2624
1.0846
λ
crb
2.0002
2.0017
0.8980
0.9028
0.6312
0.5423
15°
λ
cru
4.3946
4.4007
1.9716
1.9799
1.3859
1.1915
λ
crb
2.1154
2.1187
0.9517
0.9561
0.6683
0.5741
30°
λ
cru
5.8966
5.9317
2.6496
2.6496
1.8586
1.6020
λ
crb
2.5365
2.5491
1.1519
1.1520
0.8047
0.6909Table 7
Comparison of critical buckling loadλ
cru
=
N
x
x
cr
o
b
2
/
π
2
D
m for a simply supported FGM plate with a
/
h
= 100 and r
/
a
= 0.2. The effective material properties are computed by rule of mixtures. In order to be consistent with the literature, the properties of the metallic phase are used for normalization.
Gradient index,n
[51]
Present
% difference
0
5.2611
5.2831
−0.42
0.2
4.6564
4.6919
−0.76
1
3.6761
3.663
0.36
2
3.3672
3.3961
−0.86
5
3.1238
3.1073
0.53
10
2.9366
2.8947
1.43Figure 7
Effect of plate aspect ratioa
/
b and gradient index on the critical buckling load for a simply supported FGM plate under uniaxial compression with a
/
h
=
10.Figure 8
Variation of the critical buckling load,λ
cru
=
N
x
x
cr
o
b
2
/
π
2
D
c, with cutout dimensions for a square FGM plate with a
/
h
=
10 subjected to uniaxial compressive loading for different gradient index n and various boundary conditions.Thermal Buckling. The thermal buckling behaviour of simply supported functionally graded skew plate is studied next. The top surface is ceramic rich and the bottom surface is metal rich. The FGM plate considered here consists of aluminum and alumina. Young’s modulus, the thermal conductivity, and the coefficient of thermal expansion for alumina are E
c
= 380 GPa, κ
c
=10.4 W/mK, and α
c
=
7.4
×
10
-
6 1/°C and for aluminum E
m
= 70 GPa, κ
m
= 204 W/mK, and α
m
=
23
×
10
-
6 1/°C, respectively. Poisson’s ratio is chosen as constant, ν
= 0.3. The temperature rise of T
m
= 5°C in the metal-rich surface of the plate is assumed in the present study. In addition to nonlinear temperature distribution across the plate thickness, the linear case is also considered in the present analysis by truncating the higher order terms in (25). The plate is of uniform thickness and simply supported on all four edges. Table 8 shows the convergence of the critical buckling temperature with mesh size for different gradient index, n. It can be seen that the results from the present formulation are in very good agreement with the available solution. The influence of the plate aspect ratio a
/
b and the skew angle ψ on the critical buckling temperature for a simply supported square FGM plate are shown in Figures 9 and 10. It is seen that increasing the plate aspect ratio decreases the critical buckling temperature for both linear and nonlinear temperature distribution through the thickness. The critical buckling temperature increases with increase in the skew angle. The influence of the gradient index n is also shown in Figure 10. It is seen that with increasing gradient index, n, the critical buckling temperature decreases. This is due to the increase in the metallic volume fraction that degrades the overall stiffness of the structure. Figure 11 shows the influence of the cutout radius and the material gradient index on the critical buckling temperature. Both linear and nonlinear temperature distribution through the thickness are assumed. Again, it is seen that the combined effect of increasing the gradient index n and the cutout radius r
/
a is to lower the buckling temperature. For gradient index n
= 0, there is no difference between the linear and the nonlinear temperature distribution through the thickness as the material is homogeneous through the thickness, while, for n
>
0, the material is heterogeneous through the thickness with different thermal property.Table 8
Convergence of the critical buckling temperature for a simply supported FGM skew plate witha
/
h
= 10 and a
/
b
= 1. Nonlinear temperature rise through the thickness of the plate is assumed.
Mesh
Gradient index,n
0
1
5
10
8
×
8
3383.40
2054.61
1539.24
1496.36
16
×
16
3286.90
1995.07
1495.25
1453.99
32
×
32
3263.91
1980.96
1484.76
1443.86
40
×
40
3261.17
1979.30
1483.51
1442.60
[9]
3257.47
1977.01
1481.83
1441.02Figure 9
Critical buckling temperature as a function of plate aspect ratioa
/
b with linear and nonlinear temperature distribution through the thickness with gradient index n
= 5.Figure 10
Critical buckling temperature as a function of skew angleψ for a simply supported square FGM plate with a
/
h
=
10. Both linear and nonlinear temperature distribution through the thickness are assumed.Figure 11
Influence of cutout size on the critical buckling temperature for a square simply supported FGM plate with a centrally located circular cutout witha
/
h
= 10 for various gradient indices n. Linear and nonlinear temperature distribution through the thickness are assumed.


## 5. Conclusion

In this paper, we applied the cell-based smoothed finite element method with discrete shear gap technique to study the static and the dynamic response of functionally graded materials. The first-order shear deformation theory was used to describe the plate kinematics. The efficiency and accuracy of the present approach is demonstrated with few numerical examples. This improved finite element technique shows insensitivity to shear locking and produces excellent results in static bending, free vibration, and buckling of functionally graded plates.



---

*Source: 247932-2014-04-08.xml*