# A Novel Higher-Order Shear and Normal Deformable Plate Theory for the Static, Free Vibration and Buckling Analysis of Functionally Graded Plates

Authors: Shi-Chao Yi; Lin-Quan Yao; Bai-Jian Tang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Closed-form solution of a special higher-order shear and normal deformable plate theory is presented for the static situations, natural frequencies, and buckling responses of simple supported functionally graded materials plates (FGMs). Distinguished from the usual theories, the uniqueness is the differentia of the new plate theory. Each individual FGM plate has special characteristics, such as material properties and length-thickness ratio. These distinctive attributes determine a set of orthogonal polynomials, and then the polynomials can form an exclusive plate theory. Thus, the novel plate theory has two merits: one is the orthogonality, where the majority of the coefficients of the equations derived from Hamilton’s principle are zero; the other is the flexibility, where the order of the plate theory can be arbitrarily set. Numerical examples with different shapes of plates are presented and the achieved results are compared with the reference solutions available in the literature. Several aspects of the model involving relevant parameters, length-to-thickness, stiffness ratios, and so forth affected by static and dynamic situations are elaborate analyzed in detail. As a consequence, the applicability and the effectiveness of the present method for accurately computing deflection, stresses, natural frequencies, and buckling response of various FGM plates are demonstrated.


---

## Body

## 1. Introduction

Functionally graded materials (FGMs) have continuous transition of material properties as a function of position along certain directions and thus are regarded as most promising applications of advanced composite materials as opposed to traditional isotropic and homogeneous materials. The gradual variation of material properties can be tailored to suit specific purposes in engineering design. Design of aircraft and space vehicles structures, electronic, and biomedical installations are some examples where FGM can be fruitfully exploited. So analyzing the static, vibration, and buckling problems of this structure is particularly important [1–4].The common analysis theory is classical Kirchhoff thin plate theory (CLT) [5–7], which ignores transverse shear effects and provides reasonable results for thin plates. However, it may not obtain accurate results for moderately thick plates. A development on the CLT is the first-order shear deformation theory (FSDT) [8–11] such as the Reissner–Mindlin moderately thick plate theory which gives reasons for transverse shear effects but needs a shear correction factor. Second- and higher-order shear deformation plate theories [12–24] use higher-order polynomials to express displacement components through the plate thickness and do not require shear correction factors.And meanwhile, some novel approaches have been performed to analyze the static, dynamic, and stability behavior of functionally graded plates. Carrera unified formulation (CUF) was exploited to obtain exact Navier solutions [25] and to develop advanced finite elements [26] and hierarchical Ritz-based models [27]. Shimpi [28, 29] has developed a two variable refined plate theory, which is based on the assumption that the in-plane and transverse displacements consist of bending and shear components in which the bending components do not contribute toward shear forces and, likewise, the shear components do not contribute toward bending moments. Recently, this theory was successfully extended to FGM plate [30]. Zenkour [31] used sinusoidal shear deformation plate theory to develop a general formulation for FGM sandwich plates made of three layers with isotropic ceramic/metal FGM facings and homogeneous ceramic as core.In past two decades, more than 30 kinds of plate theories have been developed and are detailed in [32, 33]. However, choosing the suitable plate theory is still a controversial issue and the main difficulty is the diversity of the FGM plate/shell. The vast majority of researchers tend to transplant some plate theory to the FGM plates, and the gradient variety is simply treated with the integral calculation. Corresponding numerical examples verify that these kinds of approaches are feasible, but there are some matters to be resolved. First, more situations need to be considered, especially volume fraction model. That is, power model, exponential model, reciprocal mode, and so on. Second, the governing equation becomes much more complicated to cause inaccurate results. Solve these aggravating issues actually simple, only need to detail with one thing: how to transform the FGM plate into orthotropic plates or laminated plates. In fact, Batra et al. offered a hint. They used a higher-order shear and normal deformable plate theory (HOSNDPT) to analyze static deformations and free and forced vibrations of moderately thick plates [34–36]. In the HOSNDPT, three-dimensional Hooke’s law is used to derive constitutive relations for various kinetic variables in terms of the kinematic variables and the governing equation become very concise by introducing orthonormalized Legendre polynomials. For the functionally graded plates, Young’s moduli and the density distribution vary with z. Like HOSNDPT, Young’s moduli are considered as the function on z (in HOSNDPT, the function is a constant). A special set of orthogonal polynomials is determined by Young’s moduli function (also called weight function), and then these polynomials bases constitute a new higher-order shear and normal deformable plate theory. Thus, we name this type of plate theory as orthogonal higher-order shear and normal deformable plate theory (OHOSNDPT). It is obvious that the HOSNDPT with orthonormalized Legendre polynomials is a special case of OHOSNDPT.The purpose of this study is to derive the analytical solutions of orthogonal higher-order shear and normal deformable plate theory (OHOSNDPT) for bending, free vibration and buckling analysis of rectangular plates. The closed-form solutions of deflection and stress are obtained for rectangular plates with various modules by Naiver approach. The obtained results are compared with those reported in the literature. The effects of thickness ratio, modulus ratio, and load situation on deflection and stress, nature frequency, and buckling loads of FGM plates are studied.


## 2. Problem Models

Consider a rectangular plate of plan-form dimensionsa and b and uniform thickness h. The coordinate system is taken such that the x
-
y plane (z
=
0) coincides with the mid-plane of the plate (z
∈
[
-
h
/
2
,
h
/
2
]) (see Figure 1). Three different types of functionally graded plates are studied: (A) isotropic FGM plates; (B) laminated FGM plates; (C) sandwich plates with FGM skins; (D) FGM plates with arbitrary variation of stiffness; (E) laminated FGM plates.Figure 1
A sketch of the FGM plates.


### 2.1. Tape A: Isotropic FGM Plate

The plate of type A is graded from metal (bottom) to ceramic (top) (see Figure1). The volume fraction of the ceramic phase is defined as in [19](1)
V
c
=
1
2
+
z
h
p
,where z
∈
[
-
h
/
2
,
h
/
2
], h is the thickness of the plate, and p is a scalar parameter that allows the user to define gradation of material properties across the thickness direction; some references also consider this type in [25, 37–43].


### 2.2. Tape B: Sandwich Plate with FGM Core

In this type of sandwich plates, the bottom skin is isotropic (fully metal) and the top skin is isotropic (fully ceramic). The core layer is graded from metal to ceramic so that there are no interfaces between core and skins, as illustrated in Figure1. The volume fraction of the ceramic phase in the core is obtained by adapting the polynomial material law in [19](2)
V
c
=
z
c
-
h
1
h
c
p
,where z
∈
[
h
1
,
h
2
], h
c
=
h
2
-
h
1 is the thickness of the core, and p is the power-law exponent that defines the gradation of material properties across the thickness direction; many results could be referred in [37, 39, 41, 44].


### 2.3. Tape C: Sandwich Plate with FGM Skins

In C-type plates, the sandwich core is isotropic (fully ceramic) and skins are composed of a functionally graded material across the thickness direction. The bottom skin varies from a metal-rich surface (z
=
h
/
2) to a ceramic-rich surface while the top skin face varies from a ceramic-rich surface to a metal-rich surface (z
=
h
/
2), as illustrated in Figure 1. In [37, 45, 46], no interfaces need to be considered between core and skins. The volume fraction of the ceramic phase is obtained as(3)
V
c
=
z
-
h
0
h
1
-
h
0
p
z
∈
-
h
2
,
h
1
bottom  skin
,
V
c
=
1
z
∈
h
1
,
h
2
core
,
V
c
=
z
-
h
3
h
2
-
h
3
p
z
∈
h
2
,
h
2
top  skin
,where z
∈
[
-
h
/
2
,
h
/
2
] and p is a scalar parameter that allows the user to define gradation of material properties across the thickness direction of the skins.The sandwich plate C-type may be symmetric or nonsymmetric about the mid-plane as we may vary the thickness of each face. A nonsymmetric sandwich with volume fraction defined by the power-law for various exponentsp, in which top skin thickness is the same as the core thickness and the bottom skin thickness is twice the core thickness. Such thickness relation is denoted as 2-1-1. A bottom-core-top notation is being used. 1-1-1 means that skins and core have the same thickness.


### 2.4. Tape D: FGM Plates with Arbitrary Variation of Stiffness

Type-D plates is also graded from metal (bottom) to ceramic (top). But unlike Tape A (defined as power model), the volume fraction of the ceramic phase can be defined as arbitrary functions form [47]. We have the following:Power model(4)
V
c
=
z
c
-
h
1
h
c
p
,
Exponential model(5)
V
c
=
e
p
-
e
p
z
/
h
+
1
/
2
e
p
-
1
,
Reciprocal model(6)
V
c
=
2
p
+
1
2
p
1
+
2
p
z
/
h
+
1
/
2
-
1
2
p
,and so on.


### 2.5. Tape E: Laminated FGM Plates

The last tape of the plates both is considered laminated plate structure and taken the FGM of each layer into account. Then, the type-E is called laminated FGM plate. From Figure1, the volume fraction of the ceramic phase is considered as a piecewise functions and represented as(7)
V
c
=
V
c
k
z
,
p
,
z
∈
z
k
,
z
k
+
1in the kth layer, where V
c
k
(
z
,
p
) represents the kth volume fraction and V
c
k
(
z
,
p
) can be considered as any situation of Type D.For the five types of plates, A–E, the volume fraction for the metal phase is given asV
m
=
1
-
V
c. The isotropic fully ceramic plate can be seen as a particular case of plates A–E, by setting the fixed values to the parameter p in Tape A–Tape E.


## 2.1. Tape A: Isotropic FGM Plate

The plate of type A is graded from metal (bottom) to ceramic (top) (see Figure1). The volume fraction of the ceramic phase is defined as in [19](1)
V
c
=
1
2
+
z
h
p
,where z
∈
[
-
h
/
2
,
h
/
2
], h is the thickness of the plate, and p is a scalar parameter that allows the user to define gradation of material properties across the thickness direction; some references also consider this type in [25, 37–43].


## 2.2. Tape B: Sandwich Plate with FGM Core

In this type of sandwich plates, the bottom skin is isotropic (fully metal) and the top skin is isotropic (fully ceramic). The core layer is graded from metal to ceramic so that there are no interfaces between core and skins, as illustrated in Figure1. The volume fraction of the ceramic phase in the core is obtained by adapting the polynomial material law in [19](2)
V
c
=
z
c
-
h
1
h
c
p
,where z
∈
[
h
1
,
h
2
], h
c
=
h
2
-
h
1 is the thickness of the core, and p is the power-law exponent that defines the gradation of material properties across the thickness direction; many results could be referred in [37, 39, 41, 44].


## 2.3. Tape C: Sandwich Plate with FGM Skins

In C-type plates, the sandwich core is isotropic (fully ceramic) and skins are composed of a functionally graded material across the thickness direction. The bottom skin varies from a metal-rich surface (z
=
h
/
2) to a ceramic-rich surface while the top skin face varies from a ceramic-rich surface to a metal-rich surface (z
=
h
/
2), as illustrated in Figure 1. In [37, 45, 46], no interfaces need to be considered between core and skins. The volume fraction of the ceramic phase is obtained as(3)
V
c
=
z
-
h
0
h
1
-
h
0
p
z
∈
-
h
2
,
h
1
bottom  skin
,
V
c
=
1
z
∈
h
1
,
h
2
core
,
V
c
=
z
-
h
3
h
2
-
h
3
p
z
∈
h
2
,
h
2
top  skin
,where z
∈
[
-
h
/
2
,
h
/
2
] and p is a scalar parameter that allows the user to define gradation of material properties across the thickness direction of the skins.The sandwich plate C-type may be symmetric or nonsymmetric about the mid-plane as we may vary the thickness of each face. A nonsymmetric sandwich with volume fraction defined by the power-law for various exponentsp, in which top skin thickness is the same as the core thickness and the bottom skin thickness is twice the core thickness. Such thickness relation is denoted as 2-1-1. A bottom-core-top notation is being used. 1-1-1 means that skins and core have the same thickness.


## 2.4. Tape D: FGM Plates with Arbitrary Variation of Stiffness

Type-D plates is also graded from metal (bottom) to ceramic (top). But unlike Tape A (defined as power model), the volume fraction of the ceramic phase can be defined as arbitrary functions form [47]. We have the following:Power model(4)
V
c
=
z
c
-
h
1
h
c
p
,
Exponential model(5)
V
c
=
e
p
-
e
p
z
/
h
+
1
/
2
e
p
-
1
,
Reciprocal model(6)
V
c
=
2
p
+
1
2
p
1
+
2
p
z
/
h
+
1
/
2
-
1
2
p
,and so on.


## 2.5. Tape E: Laminated FGM Plates

The last tape of the plates both is considered laminated plate structure and taken the FGM of each layer into account. Then, the type-E is called laminated FGM plate. From Figure1, the volume fraction of the ceramic phase is considered as a piecewise functions and represented as(7)
V
c
=
V
c
k
z
,
p
,
z
∈
z
k
,
z
k
+
1in the kth layer, where V
c
k
(
z
,
p
) represents the kth volume fraction and V
c
k
(
z
,
p
) can be considered as any situation of Type D.For the five types of plates, A–E, the volume fraction for the metal phase is given asV
m
=
1
-
V
c. The isotropic fully ceramic plate can be seen as a particular case of plates A–E, by setting the fixed values to the parameter p in Tape A–Tape E.


## 3. Orthogonal Higher-Order Shear and Normal Deformable Plate Theory

### 3.1. Orthogonal Polynomials

Definition 1.
f
(
z
) is said to be orthogonal to g
(
z
) for the interval [
a
,
b
] with respect to the weight function ψ
(
z
) if (8)
∫
a
b
ψ
z
f
z
g
z
d
z
=
0
,and L
~
0
(
z
)
,
L
~
1
(
z
)
,
…
,
L
~
n
(
z
)
,
… is said to be an orthogonal set of functions for the interval [
a
,
b
] with respect to the weight function ψ
(
z
) if(9)
∫
a
b
ψ
z
L
~
j
z
L
~
k
z
d
z
=
0
,
j
≠
k
;
A
k
>
0
,
j
=
k
.In addition, if A
k
=
1 for each k
=
0,1
,
…
,
n
,
…, the set is said to be orthonormal.Theorem 2.
The set of polynomial functions{
L
~
i
(
z
)
}
0
∞ obtained in the following way is orthogonal on [
a
,
b
] with respect to the weight function ψ
(
z
): (10)
L
~
0
z
=
1
,
L
~
1
z
=
z
-
B
1
,
for  each
z
in
a
,
b
,where(11)
B
1
=
∫
a
b
z
ψ
z
L
~
0
z
2
d
z
∫
a
b
ψ
z
L
~
0
z
2
d
z
,and when k
>
1,(12)
L
~
k
z
=
x
-
B
k
L
~
k
-
1
z
-
C
k
L
~
k
-
2
z
,
for  each
z
in
a
,
b
,where (13)
B
k
=
∫
a
b
z
ψ
z
L
~
k
-
1
z
2
d
z
∫
a
b
ψ
z
L
~
k
-
1
z
2
d
z
,
C
k
=
∫
a
b
z
ψ
z
L
~
k
-
1
z
L
~
k
-
2
z
d
x
∫
a
b
ψ
z
L
~
k
-
2
z
2
d
z
.For arbitrarily integrable weight functionψ
(
z
), the theorem provides a recursive procedure for constructing a set of orthogonal polynomials.Remark 3.
(i) Ifψ
(
z
) is a piecewise function, the integration of ψ
(
z
) in Definition 1 and Theorem 2 can be represented as the sum of the integrations in subinterval. (ii) It is simple to show that there is a one-for-one map between the weight function ψ
(
z
) and the orthogonal polynomials {
L
~
i
(
z
)
}
0
∞. (iii) The normalized orthogonal polynomial functions set {
L
i
(
z
)
}
0
∞ are obtained by {
L
~
i
(
z
)
}
0
∞ divided by {
A
i
1
/
2
}
0
∞. That is, L
i
(
z
)
=
L
~
i
(
z
)
/
(
A
i
1
/
2
).


### 3.2. Orthogonal Higher-Order Shear and Normal Deformable Plate Theory

A 3D displacement function on the surface and thickness direction can be separated, and the thickness direction can be expanded by orthogonal polynomial{
L
i
(
z
)
}
0
∞ (Remark 3). The displacement field is assumed to be of the form(14)
u
x
,
y
,
z
,
t
=
u
x
,
y
,
z
,
t
v
x
,
y
,
z
,
t
w
x
,
y
,
z
,
t
=
∑
i
=
0
K
u
i
x
,
y
,
t
v
i
x
,
y
,
t
w
i
x
,
y
,
t
L
i
z
.When K
>
1, the plate theory is called higher-order theory.The strain-displacement relationships are given as(15)
ε
=
ε
x
ε
y
ε
z
γ
y
z
γ
z
x
γ
x
y
=
∑
i
=
0
K
∂
u
i
x
,
y
∂
x
∂
v
i
x
,
y
∂
y
∑
j
=
0
K
w
j
x
,
y
d
j
i
∂
w
i
x
,
y
∂
y
+
∑
j
=
0
K
v
j
x
,
y
d
j
i
∂
w
i
x
,
y
∂
x
+
∑
j
=
0
K
u
j
x
,
y
d
j
i
∂
v
i
x
,
y
∂
x
+
∂
u
i
x
,
y
∂
y
L
i
z
.The derivative of theith orthogonal polynomial is a polynomial of degree i
-
1, which can be linearly represented by the first i
-
1 order orthogonal polynomials. Then, it can be represented as(16)
L
i
z
=
∑
j
=
0
K
d
i
j
L
j
z
,where d
i
j is constant.


### 3.3. Constitutive Relations for a Laminates

For the functionally graded material detailed in Section2, the constitutive relations C is(17)
σ
x
x
σ
y
y
σ
z
z
τ
y
z
τ
z
x
τ
x
y
=
C
~
11
C
~
12
C
~
12
0
0
0
C
~
12
C
~
11
C
~
12
0
0
0
C
~
12
C
~
12
C
~
11
0
0
0
0
0
0
C
~
44
0
0
0
0
0
0
C
~
44
0
0
0
0
0
0
C
~
44
ε
x
x
ε
y
y
ε
z
z
γ
y
z
γ
z
x
γ
x
y
,where(18)
C
~
11
=
E
1
-
ν
2
1
-
3
ν
2
-
2
ν
3
,
C
~
12
=
E
ν
+
ν
2
1
-
3
ν
2
-
2
ν
3
,
C
~
44
=
G
,where E is the modulus of elasticity, ν is Poisson’s ratio, and G is the shear modulus G
=
E
/
2
(
1
-
ν
).Based on the volume fraction of the constituent material for Mori–Tanaka scheme, the Young modulus and density of some of FG plates can be written as a functions of thickness coordinate,z, as follows:Power model(19)
V
c
=
z
c
-
h
1
h
c
p
,
E
z
=
E
m
+
E
c
-
E
m
z
c
-
h
1
h
c
p
=
E
0
ψ
z
,
ρ
z
=
ρ
m
+
ρ
c
-
ρ
m
z
c
-
h
1
h
c
p
,
Exponential model(20)
V
c
=
e
p
-
e
p
z
/
h
+
1
/
2
e
p
-
1
,
E
z
=
E
m
e
p
z
/
h
+
1
/
2
=
E
0
ψ
z
,
ρ
z
=
ρ
m
e
p
z
/
h
+
1
/
2
,
Reciprocal model(21)
V
c
=
2
p
+
1
2
p
1
+
2
p
z
/
h
+
1
/
2
-
1
2
p
,
E
z
=
E
m
p
∗
z
/
h
+
3
/
2
=
E
0
ψ
z
,
ρ
z
=
ρ
m
p
∗
z
/
h
+
3
/
2
,where ψ
(
z
) is called stiffness function.Remark 4.
(i) The modulusE and the Poisson’s ratio ν are the functions of the z in the FGM plates, then the C
~   is related to z. If Poisson’s ratios of two different materials are the same, C
~ can be represented as C
∗
ψ
(
z
) with E
=
E
0
∗
ψ
(
z
); that is, the element of the matrix C
~, C
~
i
j, can be represented as C
i
j
∗
ψ
(
z
). (ii) ψ
(
z
) can be determined by the volume fraction of the ceramic phase in Tape A–E, and ψ
(
z
) is able to find the corresponding orthonormal polynomial functions set {
L
i
(
z
)
}
0
∞ by Theorem 2.


### 3.4. Equations of Motion

For deriving the equilibrium equations for buckling analysis using the defined displacement model, the principle of minimum potential energy (PMPE) is opted due to its simplicity and also because its application gives simultaneously the natural boundary conditions that are to be used with theory. In analytical form, it can be written in [50](22)
0
=
δ
U
+
K
+
∫
Ω
0
∫
-
h
/
2
h
/
2
σ
^
x
x
δ
ε
^
x
x
+
σ
^
y
y
δ
ε
^
y
y
+
σ
^
x
y
δ
ε
^
x
y
d
z
d
x
d
y
,where U is the total strain energy due to deformation, the integral part in (22) is the potential of the external loads, the second term in the above equation is the potential energy due to the in-plane stresses σ
^
x
x, σ
^
y
y, and σ
^
x
y produced due to applied middle plane loads and ε
^
x
x, ε
^
y
y, and ε
^
x
y are in-plane strains produced by transverse deflection w and δ denotes the variational symbol. Substituting the appropriate energy expression in the above equation, the final expression can thus be written as(23)
0
=
∫
Ω
0
∫
-
h
/
2
h
/
2
σ
x
x
δ
ε
x
x
+
σ
y
y
δ
ε
y
y
+
σ
z
z
δ
ε
z
z
+
σ
y
z
δ
ε
y
z
+
σ
z
x
δ
ε
z
x
+
σ
x
y
δ
ε
x
y
+
σ
^
x
x
δ
ε
^
x
x
+
σ
^
y
y
δ
ε
^
y
y
+
σ
^
x
y
δ
ε
^
x
y
+
ρ
u
˙
δ
u
˙
+
v
˙
δ
v
˙
+
w
˙
δ
w
˙
d
z
d
x
d
y
+
∫
Ω
0
q
b
x
,
y
δ
w
x
,
y
,
h
2
+
q
t
x
,
y
δ
w
x
,
y
,
-
h
2
d
x
d
y
,where q
b is the distributed force at the bottom (z
=
h
/
2) of the laminate, q
t is the distributed force at the top (z
=
-
h
/
2) of the laminate as shown in Figure 2, (σ
^
x
x
,
σ
^
y
y
,
σ
^
x
y
,
σ
^
y
x) are the specified stress components on the portion Γ
σ. Using (15), (23), and integrating the resulting expression by parts, and collecting the coefficients δ
u
0
,
δ
v
0
,
δ
w
0
,
…
,
δ
u
i, δ
v
i
,
δ
w
i
,
…
,
δ
u
K
,
δ
v
K
,
δ
w
K.Figure 2
Plate and coordinate system.For anyi, (
δ
u
i
,
δ
v
i
,
δ
w
i
), the following equations of equilibrium are obtained:(24)
δ
u
i
:
N
x
x
,
x
i
+
N
x
y
,
y
i
+
Q
z
x
i
=
∑
j
=
0
K
I
i
j
u
¨
j
,
δ
v
i
:
N
x
y
,
x
i
+
N
y
y
,
y
i
+
Q
y
z
i
=
∑
j
=
0
K
I
i
j
v
¨
j
,
δ
w
i
:
N
z
x
,
x
i
+
N
y
z
,
y
i
+
Q
z
z
i
+
∑
j
=
0
K
N
^
x
x
i
j
∂
2
w
j
∂
x
2
+
∑
j
=
0
K
N
^
y
y
i
j
∂
2
w
j
∂
y
2
+
∑
j
=
0
K
N
^
x
y
i
j
∂
2
w
j
∂
x
∂
y
+
q
b
x
,
y
L
i
h
2
+
q
t
x
,
y
L
i
-
h
2
=
∑
j
=
0
K
I
i
j
w
¨
j
,where (25)
N
12
i
j
=
∫
-
h
/
2
h
/
2
σ
12
L
i
z
L
j
z
d
z
,
N
^
12
i
j
=
∫
-
h
/
2
h
/
2
σ
^
12
L
i
z
L
j
z
d
z
,
1,2
=
x
,
y
,
z
,
Q
i
j
i
=
∫
-
h
/
2
h
/
2
∂
L
i
z
∂
z
σ
i
j
d
z
,
I
i
j
=
∫
-
h
/
2
h
/
2
ρ
L
i
z
L
j
z
d
z
.Substituting equation (15) into equation (24) with stress-strain relations equation (17) and equation (18), the governing equations of motion on (
δ
u
i
,
δ
v
i
,
δ
w
i
) are obtained as(26)
δ
u
i
:
C
11
∂
2
u
i
x
,
y
∂
x
2
+
C
12
∂
2
v
i
x
,
y
∂
x
∂
y
+
C
13
∑
j
=
0
K
d
j
i
∂
w
j
x
,
y
∂
x
+
C
66
∂
∂
v
i
x
,
y
/
∂
x
+
∂
u
i
x
,
y
/
∂
y
∂
y
-
C
55
∑
j
=
0
K
d
i
j
∂
w
j
x
,
y
∂
x
+
∑
k
=
0
K
d
i
k
d
j
k
u
j
x
,
y
=
∑
j
=
0
K
I
i
j
u
¨
j
,
δ
v
i
:
C
21
∂
2
u
i
x
,
y
∂
x
∂
y
+
C
22
∂
2
v
i
x
,
y
∂
y
2
+
C
23
∑
j
=
0
K
d
j
i
∂
w
j
x
,
y
∂
y
+
C
66
∂
∂
v
i
x
,
y
/
∂
x
+
∂
u
i
x
,
y
/
∂
y
∂
x
-
C
44
∑
j
=
0
K
d
i
j
∂
w
j
x
,
y
∂
y
+
∑
k
=
0
K
d
i
k
d
j
k
v
j
x
,
y
=
∑
j
=
0
K
I
i
j
v
¨
j
,
δ
w
i
:
C
55
∂
2
w
i
x
,
y
∂
x
2
+
C
55
∑
j
=
0
K
d
j
i
∂
u
i
x
,
y
∂
x
+
C
44
∂
2
w
i
x
,
y
∂
y
2
+
C
44
∑
j
=
0
K
d
j
i
∂
v
i
x
,
y
∂
y
+
C
31
∑
j
=
0
K
d
i
j
∂
u
j
x
,
y
∂
x
+
C
32
∑
j
=
0
K
d
i
j
∂
v
j
x
,
y
∂
y
+
C
33
∑
j
=
0
K
∑
k
=
0
K
d
i
k
d
j
k
w
j
x
,
y
+
∑
j
=
0
K
N
^
x
x
i
j
∂
2
w
j
∂
x
2
+
∑
j
=
0
K
N
^
y
y
i
j
∂
2
w
j
∂
y
2
+
∑
j
=
0
K
N
^
x
y
i
j
∂
2
w
j
∂
x
∂
y
+
q
b
x
,
y
L
i
h
2
+
q
t
x
,
y
L
i
-
h
2
=
∑
j
=
0
K
I
i
j
w
¨
j
.


## 3.1. Orthogonal Polynomials

Definition 1.
f
(
z
) is said to be orthogonal to g
(
z
) for the interval [
a
,
b
] with respect to the weight function ψ
(
z
) if (8)
∫
a
b
ψ
z
f
z
g
z
d
z
=
0
,and L
~
0
(
z
)
,
L
~
1
(
z
)
,
…
,
L
~
n
(
z
)
,
… is said to be an orthogonal set of functions for the interval [
a
,
b
] with respect to the weight function ψ
(
z
) if(9)
∫
a
b
ψ
z
L
~
j
z
L
~
k
z
d
z
=
0
,
j
≠
k
;
A
k
>
0
,
j
=
k
.In addition, if A
k
=
1 for each k
=
0,1
,
…
,
n
,
…, the set is said to be orthonormal.Theorem 2.
The set of polynomial functions{
L
~
i
(
z
)
}
0
∞ obtained in the following way is orthogonal on [
a
,
b
] with respect to the weight function ψ
(
z
): (10)
L
~
0
z
=
1
,
L
~
1
z
=
z
-
B
1
,
for  each
z
in
a
,
b
,where(11)
B
1
=
∫
a
b
z
ψ
z
L
~
0
z
2
d
z
∫
a
b
ψ
z
L
~
0
z
2
d
z
,and when k
>
1,(12)
L
~
k
z
=
x
-
B
k
L
~
k
-
1
z
-
C
k
L
~
k
-
2
z
,
for  each
z
in
a
,
b
,where (13)
B
k
=
∫
a
b
z
ψ
z
L
~
k
-
1
z
2
d
z
∫
a
b
ψ
z
L
~
k
-
1
z
2
d
z
,
C
k
=
∫
a
b
z
ψ
z
L
~
k
-
1
z
L
~
k
-
2
z
d
x
∫
a
b
ψ
z
L
~
k
-
2
z
2
d
z
.For arbitrarily integrable weight functionψ
(
z
), the theorem provides a recursive procedure for constructing a set of orthogonal polynomials.Remark 3.
(i) Ifψ
(
z
) is a piecewise function, the integration of ψ
(
z
) in Definition 1 and Theorem 2 can be represented as the sum of the integrations in subinterval. (ii) It is simple to show that there is a one-for-one map between the weight function ψ
(
z
) and the orthogonal polynomials {
L
~
i
(
z
)
}
0
∞. (iii) The normalized orthogonal polynomial functions set {
L
i
(
z
)
}
0
∞ are obtained by {
L
~
i
(
z
)
}
0
∞ divided by {
A
i
1
/
2
}
0
∞. That is, L
i
(
z
)
=
L
~
i
(
z
)
/
(
A
i
1
/
2
).


## 3.2. Orthogonal Higher-Order Shear and Normal Deformable Plate Theory

A 3D displacement function on the surface and thickness direction can be separated, and the thickness direction can be expanded by orthogonal polynomial{
L
i
(
z
)
}
0
∞ (Remark 3). The displacement field is assumed to be of the form(14)
u
x
,
y
,
z
,
t
=
u
x
,
y
,
z
,
t
v
x
,
y
,
z
,
t
w
x
,
y
,
z
,
t
=
∑
i
=
0
K
u
i
x
,
y
,
t
v
i
x
,
y
,
t
w
i
x
,
y
,
t
L
i
z
.When K
>
1, the plate theory is called higher-order theory.The strain-displacement relationships are given as(15)
ε
=
ε
x
ε
y
ε
z
γ
y
z
γ
z
x
γ
x
y
=
∑
i
=
0
K
∂
u
i
x
,
y
∂
x
∂
v
i
x
,
y
∂
y
∑
j
=
0
K
w
j
x
,
y
d
j
i
∂
w
i
x
,
y
∂
y
+
∑
j
=
0
K
v
j
x
,
y
d
j
i
∂
w
i
x
,
y
∂
x
+
∑
j
=
0
K
u
j
x
,
y
d
j
i
∂
v
i
x
,
y
∂
x
+
∂
u
i
x
,
y
∂
y
L
i
z
.The derivative of theith orthogonal polynomial is a polynomial of degree i
-
1, which can be linearly represented by the first i
-
1 order orthogonal polynomials. Then, it can be represented as(16)
L
i
z
=
∑
j
=
0
K
d
i
j
L
j
z
,where d
i
j is constant.


## 3.3. Constitutive Relations for a Laminates

For the functionally graded material detailed in Section2, the constitutive relations C is(17)
σ
x
x
σ
y
y
σ
z
z
τ
y
z
τ
z
x
τ
x
y
=
C
~
11
C
~
12
C
~
12
0
0
0
C
~
12
C
~
11
C
~
12
0
0
0
C
~
12
C
~
12
C
~
11
0
0
0
0
0
0
C
~
44
0
0
0
0
0
0
C
~
44
0
0
0
0
0
0
C
~
44
ε
x
x
ε
y
y
ε
z
z
γ
y
z
γ
z
x
γ
x
y
,where(18)
C
~
11
=
E
1
-
ν
2
1
-
3
ν
2
-
2
ν
3
,
C
~
12
=
E
ν
+
ν
2
1
-
3
ν
2
-
2
ν
3
,
C
~
44
=
G
,where E is the modulus of elasticity, ν is Poisson’s ratio, and G is the shear modulus G
=
E
/
2
(
1
-
ν
).Based on the volume fraction of the constituent material for Mori–Tanaka scheme, the Young modulus and density of some of FG plates can be written as a functions of thickness coordinate,z, as follows:Power model(19)
V
c
=
z
c
-
h
1
h
c
p
,
E
z
=
E
m
+
E
c
-
E
m
z
c
-
h
1
h
c
p
=
E
0
ψ
z
,
ρ
z
=
ρ
m
+
ρ
c
-
ρ
m
z
c
-
h
1
h
c
p
,
Exponential model(20)
V
c
=
e
p
-
e
p
z
/
h
+
1
/
2
e
p
-
1
,
E
z
=
E
m
e
p
z
/
h
+
1
/
2
=
E
0
ψ
z
,
ρ
z
=
ρ
m
e
p
z
/
h
+
1
/
2
,
Reciprocal model(21)
V
c
=
2
p
+
1
2
p
1
+
2
p
z
/
h
+
1
/
2
-
1
2
p
,
E
z
=
E
m
p
∗
z
/
h
+
3
/
2
=
E
0
ψ
z
,
ρ
z
=
ρ
m
p
∗
z
/
h
+
3
/
2
,where ψ
(
z
) is called stiffness function.Remark 4.
(i) The modulusE and the Poisson’s ratio ν are the functions of the z in the FGM plates, then the C
~   is related to z. If Poisson’s ratios of two different materials are the same, C
~ can be represented as C
∗
ψ
(
z
) with E
=
E
0
∗
ψ
(
z
); that is, the element of the matrix C
~, C
~
i
j, can be represented as C
i
j
∗
ψ
(
z
). (ii) ψ
(
z
) can be determined by the volume fraction of the ceramic phase in Tape A–E, and ψ
(
z
) is able to find the corresponding orthonormal polynomial functions set {
L
i
(
z
)
}
0
∞ by Theorem 2.


## 3.4. Equations of Motion

For deriving the equilibrium equations for buckling analysis using the defined displacement model, the principle of minimum potential energy (PMPE) is opted due to its simplicity and also because its application gives simultaneously the natural boundary conditions that are to be used with theory. In analytical form, it can be written in [50](22)
0
=
δ
U
+
K
+
∫
Ω
0
∫
-
h
/
2
h
/
2
σ
^
x
x
δ
ε
^
x
x
+
σ
^
y
y
δ
ε
^
y
y
+
σ
^
x
y
δ
ε
^
x
y
d
z
d
x
d
y
,where U is the total strain energy due to deformation, the integral part in (22) is the potential of the external loads, the second term in the above equation is the potential energy due to the in-plane stresses σ
^
x
x, σ
^
y
y, and σ
^
x
y produced due to applied middle plane loads and ε
^
x
x, ε
^
y
y, and ε
^
x
y are in-plane strains produced by transverse deflection w and δ denotes the variational symbol. Substituting the appropriate energy expression in the above equation, the final expression can thus be written as(23)
0
=
∫
Ω
0
∫
-
h
/
2
h
/
2
σ
x
x
δ
ε
x
x
+
σ
y
y
δ
ε
y
y
+
σ
z
z
δ
ε
z
z
+
σ
y
z
δ
ε
y
z
+
σ
z
x
δ
ε
z
x
+
σ
x
y
δ
ε
x
y
+
σ
^
x
x
δ
ε
^
x
x
+
σ
^
y
y
δ
ε
^
y
y
+
σ
^
x
y
δ
ε
^
x
y
+
ρ
u
˙
δ
u
˙
+
v
˙
δ
v
˙
+
w
˙
δ
w
˙
d
z
d
x
d
y
+
∫
Ω
0
q
b
x
,
y
δ
w
x
,
y
,
h
2
+
q
t
x
,
y
δ
w
x
,
y
,
-
h
2
d
x
d
y
,where q
b is the distributed force at the bottom (z
=
h
/
2) of the laminate, q
t is the distributed force at the top (z
=
-
h
/
2) of the laminate as shown in Figure 2, (σ
^
x
x
,
σ
^
y
y
,
σ
^
x
y
,
σ
^
y
x) are the specified stress components on the portion Γ
σ. Using (15), (23), and integrating the resulting expression by parts, and collecting the coefficients δ
u
0
,
δ
v
0
,
δ
w
0
,
…
,
δ
u
i, δ
v
i
,
δ
w
i
,
…
,
δ
u
K
,
δ
v
K
,
δ
w
K.Figure 2
Plate and coordinate system.For anyi, (
δ
u
i
,
δ
v
i
,
δ
w
i
), the following equations of equilibrium are obtained:(24)
δ
u
i
:
N
x
x
,
x
i
+
N
x
y
,
y
i
+
Q
z
x
i
=
∑
j
=
0
K
I
i
j
u
¨
j
,
δ
v
i
:
N
x
y
,
x
i
+
N
y
y
,
y
i
+
Q
y
z
i
=
∑
j
=
0
K
I
i
j
v
¨
j
,
δ
w
i
:
N
z
x
,
x
i
+
N
y
z
,
y
i
+
Q
z
z
i
+
∑
j
=
0
K
N
^
x
x
i
j
∂
2
w
j
∂
x
2
+
∑
j
=
0
K
N
^
y
y
i
j
∂
2
w
j
∂
y
2
+
∑
j
=
0
K
N
^
x
y
i
j
∂
2
w
j
∂
x
∂
y
+
q
b
x
,
y
L
i
h
2
+
q
t
x
,
y
L
i
-
h
2
=
∑
j
=
0
K
I
i
j
w
¨
j
,where (25)
N
12
i
j
=
∫
-
h
/
2
h
/
2
σ
12
L
i
z
L
j
z
d
z
,
N
^
12
i
j
=
∫
-
h
/
2
h
/
2
σ
^
12
L
i
z
L
j
z
d
z
,
1,2
=
x
,
y
,
z
,
Q
i
j
i
=
∫
-
h
/
2
h
/
2
∂
L
i
z
∂
z
σ
i
j
d
z
,
I
i
j
=
∫
-
h
/
2
h
/
2
ρ
L
i
z
L
j
z
d
z
.Substituting equation (15) into equation (24) with stress-strain relations equation (17) and equation (18), the governing equations of motion on (
δ
u
i
,
δ
v
i
,
δ
w
i
) are obtained as(26)
δ
u
i
:
C
11
∂
2
u
i
x
,
y
∂
x
2
+
C
12
∂
2
v
i
x
,
y
∂
x
∂
y
+
C
13
∑
j
=
0
K
d
j
i
∂
w
j
x
,
y
∂
x
+
C
66
∂
∂
v
i
x
,
y
/
∂
x
+
∂
u
i
x
,
y
/
∂
y
∂
y
-
C
55
∑
j
=
0
K
d
i
j
∂
w
j
x
,
y
∂
x
+
∑
k
=
0
K
d
i
k
d
j
k
u
j
x
,
y
=
∑
j
=
0
K
I
i
j
u
¨
j
,
δ
v
i
:
C
21
∂
2
u
i
x
,
y
∂
x
∂
y
+
C
22
∂
2
v
i
x
,
y
∂
y
2
+
C
23
∑
j
=
0
K
d
j
i
∂
w
j
x
,
y
∂
y
+
C
66
∂
∂
v
i
x
,
y
/
∂
x
+
∂
u
i
x
,
y
/
∂
y
∂
x
-
C
44
∑
j
=
0
K
d
i
j
∂
w
j
x
,
y
∂
y
+
∑
k
=
0
K
d
i
k
d
j
k
v
j
x
,
y
=
∑
j
=
0
K
I
i
j
v
¨
j
,
δ
w
i
:
C
55
∂
2
w
i
x
,
y
∂
x
2
+
C
55
∑
j
=
0
K
d
j
i
∂
u
i
x
,
y
∂
x
+
C
44
∂
2
w
i
x
,
y
∂
y
2
+
C
44
∑
j
=
0
K
d
j
i
∂
v
i
x
,
y
∂
y
+
C
31
∑
j
=
0
K
d
i
j
∂
u
j
x
,
y
∂
x
+
C
32
∑
j
=
0
K
d
i
j
∂
v
j
x
,
y
∂
y
+
C
33
∑
j
=
0
K
∑
k
=
0
K
d
i
k
d
j
k
w
j
x
,
y
+
∑
j
=
0
K
N
^
x
x
i
j
∂
2
w
j
∂
x
2
+
∑
j
=
0
K
N
^
y
y
i
j
∂
2
w
j
∂
y
2
+
∑
j
=
0
K
N
^
x
y
i
j
∂
2
w
j
∂
x
∂
y
+
q
b
x
,
y
L
i
h
2
+
q
t
x
,
y
L
i
-
h
2
=
∑
j
=
0
K
I
i
j
w
¨
j
.


## 4. Analytical Solution

See Figure3; the rectangular plate with length a and width b under consideration is solved for the following simply supported boundary conditions prescribed at all four edges:(27)
u
i
x
,
0
=
u
i
x
,
b
=
v
i
0
,
y
=
v
i
a
,
y
=
0
,
w
i
x
,
0
=
w
i
x
,
b
=
w
i
0
,
y
=
w
i
a
,
y
=
0
,
N
y
y
i
x
,
0
=
N
y
y
i
x
,
b
=
N
x
x
i
0
,
y
=
N
x
x
i
a
,
y
=
0
.Figure 3
Simply supported boundary conditions for a rectangular plate.For the analytical solution of the partial differential equations (26), the Navier method, based on double Fourier series, is used under the specified boundary conditions (27). Using Navier’s procedure, the solution of the displacement variables satisfying the above boundary conditions can be expressed in the following Fourier series:(28)
u
i
x
,
y
,
t
=
∑
m
=
1
∞
∑
n
=
1
∞
U
i
,
m
n
cos
⁡
α
x
sin
⁡
β
y
e
i
ω
t
,
0
≤
x
≤
a
;
0
≤
y
≤
b
;
i
=
0
,
…
,
K
,
v
i
x
,
y
,
t
=
∑
m
=
1
∞
∑
n
=
1
∞
V
i
,
m
n
sin
⁡
α
x
cos
⁡
β
y
e
i
ω
t
,
0
≤
x
≤
a
;
0
≤
y
≤
b
;
i
=
0
,
…
,
K
,
w
i
x
,
y
,
t
=
∑
m
=
1
∞
∑
n
=
1
∞
W
i
,
m
n
sin
⁡
α
x
sin
⁡
β
y
e
i
ω
t
,
0
≤
x
≤
a
;
0
≤
y
≤
b
;
i
=
0
,
…
,
K
,where (29)
α
=
m
π
a
,
β
=
n
π
b
,and U
i
,
m
n, V
i
,
m
n, and W
i
,
m
n are arbitrary parameters to be determined and ω is the natural frequency.


### 4.1. Static Problem

Substituting equations (28) into equations (26), the following equations are obtained for any fixed value of m and n:(30)
K
U
m
n
=
F
m
n
,where [
K
]
=
[
K
a
]
+
[
K
b
] and [
K
a
] and [
K
b
] consists of 3
×
3 matrices [
K
i
j
a
] and [
K
i
j
b
], respectively. Note that these matrices are all diagonal block matrices. [
K
i
i
a
] and [
K
i
j
b
] are represented as(31)
K
i
i
a
=
C
11
α
2
+
C
66
β
2
C
12
+
C
66
α
β
0
C
21
α
β
+
C
66
α
β
C
22
β
2
+
C
66
α
2
0
0
0
2
C
45
α
β
+
C
55
α
2
+
C
44
β
2
,
K
i
j
b
=
-
C
55
∑
k
=
0
K
d
i
k
d
j
k
-
C
45
∑
k
=
0
K
d
i
k
d
j
k
C
13
d
j
i
α
-
C
45
d
i
j
β
-
C
55
d
i
j
α
-
C
45
∑
k
=
0
K
d
i
k
d
j
k
-
C
44
∑
k
=
0
K
d
i
k
d
j
k
C
23
d
j
i
β
-
C
44
d
i
j
β
-
C
45
d
i
j
α
C
55
d
j
i
α
+
C
45
d
j
i
β
+
C
31
d
i
j
α
C
45
d
j
i
α
+
C
44
d
j
i
α
+
C
32
d
i
j
β
C
33
∑
k
=
0
K
d
i
k
d
j
k
,
U
m
n
=
U
1
,
m
n
,
V
1
,
m
n
,
W
1
,
m
n
,
U
i
,
m
n
,
V
i
,
m
n
,
W
i
,
m
n
,
U
K
,
m
n
,
V
K
,
m
n
,
W
K
,
m
n
T
,
F
=
0,0
,
F
1
,
m
n
,
0,0
,
F
i
,
m
n
,
0,0
,
F
K
,
m
n
T
,
F
i
,
m
n
=
4
a
b
∫
0
a
∫
0
b
q
b
x
,
y
L
i
h
2
+
q
t
x
,
y
L
i
-
h
2
sin
⁡
α
x
sin
⁡
β
y
d
x
d
y
.


### 4.2. Buckling Problem

Substituting equation (28) into equations of motion equation (26) and only considering N
x
x
00 and N
y
y
00 are nonzero, we get below equation for any fixed value of m and n, for buckling vibration problem:(32)
K
~
U
m
n
=
0
,where [
K
~
]
=
[
K
~
a
]
+
[
K
b
]. The [
K
~
i
i
a
] is represented as(33)
K
i
i
a
=
C
11
α
2
+
C
66
β
2
C
12
+
C
66
α
β
0
C
21
α
β
+
C
66
α
β
C
22
β
2
+
C
66
α
2
0
0
0
C
55
α
2
+
C
44
β
2
+
N
x
x
00
α
2
+
N
y
y
00
β
2
δ
i
0
.The solutions of N
x
x
00 and N
y
y
00 are obtained by solving the determinant of K
~ equal to zero.


### 4.3. Free Vibrations Problems

Substituting equation (28) into equations of motion equation (26) neglecting N
x
x
00, N
y
y
00, and external force F
m
n, we get below eigenvalue equation for any fixed value of m and n, for free vibration problem:(34)
K
-
ω
2
M
U
m
n
=
0
,where [
M
] is also a block matrix and the [
M
] is made up of 3
×
3 matrices [
M
i
j
], where M
i
j
=
I
j
i
E and E is an 3
×
3 identity matrix.


## 4.1. Static Problem

Substituting equations (28) into equations (26), the following equations are obtained for any fixed value of m and n:(30)
K
U
m
n
=
F
m
n
,where [
K
]
=
[
K
a
]
+
[
K
b
] and [
K
a
] and [
K
b
] consists of 3
×
3 matrices [
K
i
j
a
] and [
K
i
j
b
], respectively. Note that these matrices are all diagonal block matrices. [
K
i
i
a
] and [
K
i
j
b
] are represented as(31)
K
i
i
a
=
C
11
α
2
+
C
66
β
2
C
12
+
C
66
α
β
0
C
21
α
β
+
C
66
α
β
C
22
β
2
+
C
66
α
2
0
0
0
2
C
45
α
β
+
C
55
α
2
+
C
44
β
2
,
K
i
j
b
=
-
C
55
∑
k
=
0
K
d
i
k
d
j
k
-
C
45
∑
k
=
0
K
d
i
k
d
j
k
C
13
d
j
i
α
-
C
45
d
i
j
β
-
C
55
d
i
j
α
-
C
45
∑
k
=
0
K
d
i
k
d
j
k
-
C
44
∑
k
=
0
K
d
i
k
d
j
k
C
23
d
j
i
β
-
C
44
d
i
j
β
-
C
45
d
i
j
α
C
55
d
j
i
α
+
C
45
d
j
i
β
+
C
31
d
i
j
α
C
45
d
j
i
α
+
C
44
d
j
i
α
+
C
32
d
i
j
β
C
33
∑
k
=
0
K
d
i
k
d
j
k
,
U
m
n
=
U
1
,
m
n
,
V
1
,
m
n
,
W
1
,
m
n
,
U
i
,
m
n
,
V
i
,
m
n
,
W
i
,
m
n
,
U
K
,
m
n
,
V
K
,
m
n
,
W
K
,
m
n
T
,
F
=
0,0
,
F
1
,
m
n
,
0,0
,
F
i
,
m
n
,
0,0
,
F
K
,
m
n
T
,
F
i
,
m
n
=
4
a
b
∫
0
a
∫
0
b
q
b
x
,
y
L
i
h
2
+
q
t
x
,
y
L
i
-
h
2
sin
⁡
α
x
sin
⁡
β
y
d
x
d
y
.


## 4.2. Buckling Problem

Substituting equation (28) into equations of motion equation (26) and only considering N
x
x
00 and N
y
y
00 are nonzero, we get below equation for any fixed value of m and n, for buckling vibration problem:(32)
K
~
U
m
n
=
0
,where [
K
~
]
=
[
K
~
a
]
+
[
K
b
]. The [
K
~
i
i
a
] is represented as(33)
K
i
i
a
=
C
11
α
2
+
C
66
β
2
C
12
+
C
66
α
β
0
C
21
α
β
+
C
66
α
β
C
22
β
2
+
C
66
α
2
0
0
0
C
55
α
2
+
C
44
β
2
+
N
x
x
00
α
2
+
N
y
y
00
β
2
δ
i
0
.The solutions of N
x
x
00 and N
y
y
00 are obtained by solving the determinant of K
~ equal to zero.


## 4.3. Free Vibrations Problems

Substituting equation (28) into equations of motion equation (26) neglecting N
x
x
00, N
y
y
00, and external force F
m
n, we get below eigenvalue equation for any fixed value of m and n, for free vibration problem:(34)
K
-
ω
2
M
U
m
n
=
0
,where [
M
] is also a block matrix and the [
M
] is made up of 3
×
3 matrices [
M
i
j
], where M
i
j
=
I
j
i
E and E is an 3
×
3 identity matrix.


## 5. Numerical Examples

In this section, various numerical examples are presented and discussed to verify the accuracy of present solutions. Comparison studies are carried out for a large number of plates with different values of aspect ratio, thickness ratio, and various combinations of boundary conditions. The following material properties are used.Poisson’s ratio is considered constant for both materialsν
m
=
ν
c
=
ν
=
0.3.


### 5.1. Bending

In the following static examples, we consider that the plate is subjected to a bisinusoidal transverse mechanical load of amplitude loadq
t
x
,
y
=
q
¯
t
sin
⁡
π
x
/
a
sin
⁡
π
y
/
b applied at the top of the plate with q
¯
t
=
1. It should be noted that the load is applied at the top surface (
z
=
h
/
2
) as shown in Figure 2.Due to the symmetry of the problem, the displacementv, stresses σ
y, and transverse shear stresses τ
z
y are not described since their distributions along the plate thickness direction are similar to those of u
,
σ
x and τ
z
x, respectively. Considering σ
x is the first-order derivative function of the displacement u, only the stress σ
x should be discussed in this paper.


#### 5.1.1. Isotropic FGM Square Plate

In this example, an isotropic FGM square plate of type-A is considered. The plate is graded from aluminumE
m
=
70
G
P
a at the bottom to alumina E
c
=
380
G
P
a at the top. The transverse displacement, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as(35)
w
¯
=
10
h
3
E
c
q
a
4
w
,
σ
¯
x
x
=
h
2
q
a
2
σ
x
x
.An initial convergence study was performed forσ
x
x
h
/
3 and transverse displacement w
(
0
) at the center of the plate, considering p
=
1, a
/
h
=
10, and the order K of the plate theory of 1, 3, 5, and 7. Compared with the references [37] and [39], results are presented in Tables 1 and 2. As seen in these tables, it is sufficient to use K
=
5 (we consider that K
=
5 should be used in the following computation in this section).Table 1
w
- convergence study for the bending analysis of A-type plate using higher-order plate theory, p
=
1, and a
/
h
=
10.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT(
k
=
5
/
6
) [37]
[39]
Meshfree [37]
w
-
0.5846
0.5891
0.5891
0.5891
0.5889
0.5890
0.5868Table 2
σ
-
x
x convergence study for the bending analysis of A-type plate using higher-order plate theory, p
=
1, and a
/
h
=
10.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (k
=
5
/
6) [37]
[39]
Meshfree [37]
σ
-
x
x
1.4968
1.4904
1.4908
1.4908
2.0150
1.5064
1.4917In Table3, we present results for σ
¯
x
x and transverse displacement for various exponents p of the power-law (18) considering the order of the plate theory equal 7. The considered side-to-thickness ratios (
a
/
h
) are 4, 10, and 100, which means thickness h equals 0.25, 0.1, and 0.01, respectively. Results are compared with the classical plate theory (CLPT), the first-order shear deformation theory (FSDT) with a correction factor k
=
5
/
6, and those from Zenkour’s generalized shear deformation theory [38], and those from Carrera et al. [25, 39], and Neves et al. [40].Table 3
A-type plate in bending. Theσ
z
x and w according to present higher-order plate theory.
p
a
/
h
σ
x
x
(
h
/
3
)
w
(
0
)
4
10
100
4
10
100
0
Meshfree [37]
≠0
0.5278
1.3176
13.161
0.3665
0.2942
0.2803
Present
0.5144
1.3123
13.1713
0.3800
0.2961
0.2804
0.5
Meshfree [37]
≠0
0.5860
1.4680
14.673
0.5493
0.4548
0.4365
Present
≠0
0.5724
1.4591
14.6446
0.5678
0.4538
0.4326
1
Ref. [41]
≠0
0.6221
14.969
14.969
0.7171
0.5875
0.5625
CLPT [37]
0
0.8060
2.0150
20.150
0.5623
0.5623
0.5623
FSDT (k
=
5
/
6) [37]
0
0.8060
2.0150
20.150
0.7291
0.5889
0.5625
Ref. [39] N
=
4
≠0
0.6221
1.5064
14.969
0.7171
0.5875
0.5625
Ref. [40]
≠0
0.5925
1.4945
14.969
0.6997
0.5845
0.5624
Meshfree [37]
≠0
0.5911
1.4917
14.945
0.7020
0.5868
0.5647
Present
0.5841
1.4908
14.9677
0.7318
0.5891
0.5625
4
Ref. [41]
≠0
0.4877
1.1971
11.923
1.1585
0.8821
0.8286
CLPT [37]
0
0.6420
1.6049
16.049
0.8281
0.8281
0.8281
FSDT (k
=
5
/
6) [37]
0
0.6420
1.6049
16.049
1.1125
0.8736
0.828
Ref. [39] N
=
4
≠0
0.4877
1.1971
11.923
1.1585
0.8821
0.8286
Ref. [40]
≠0
0.2723
0.2778
0.2785
1.0391
0.8202
0.7784
Meshfree [37]
≠0
0.4330
1.1588
11.737
1.1108
0.8700
0.8240
Present
0.4451
1.1787
11.9209
1.1732
0.8831
0.8287
10
Ref. [41]
≠0
0.3695
0.8965
8.9077
1.3745
1.0072
0.9361
CLPT [37]
0
0.4796
1.1990
11.990
0.9354
0.9354
0.9354
FSDT (k
=
5
/
6) [37]
0
0.4796
1.1990
11.990
1.3178
0.9966
0.9360
Ref. [39] N
=
4
≠0
0.1478
0.8965
8.9077
1.3745
1.0072
0.9361
Ref. [40]
≠0
0.3227
1.1783
11.932
1.3490
0.8750
0.8286
Meshfree [37]
≠0
0.3097
0.8462
8.6010
1.3334
0.9888
0.9227
Present
≠0
0.3238
0.8766
8.9058
1.4026
1.0094
0.9362Date shown in Table3 clearly illustrated, a subtle difference appears in the moderately thick plate between them and the present higher-order plate theory are in good agreement with those from the Refs. When the side-to-thickness ratio is 100, the values, whether it is on transverse displacements or stresses, are almost the same as the values of [50].In Figure4, we present the evolution of the displacement and stresses across the thickness direction according to present shear deformation theory for various values of the exponent p and side-to-thickness ratio a
/
h
=
4, using the 5-order plate theory. Compared with the results of the meshfree method in [44], it shows that the values of the in-plane stress σ
x
x, σ
x
y and out-plane stress σ
z
x are consistent with the reference value. Strangeness is the presence of a significant difference between the present theory (K
=
5) and the meshless solutions on the longitudinal displacement w on p
=
10. From Figure 4, we find out the relative error of displacement w can up to 9% and the maximum relative error of stresses between them is less than 3%. This discord phenomenon is due to the computational error and the instability of the numerical method. On the other hand, it demonstrates the analytical solution of the orthogonal higher-order shear deformable plate theory is necessary and feasible.Figure 4
A-type square plate subjected to sinusoidal load at the top, witha
/
h
=
4, nondimensional displacement (w
¯) and stresses through the thickness to present OHOSNDPT for different values of p.In summary, it can be concluded that the present higher-order theory provides excellent solution for the isotropic FGM plates.


#### 5.1.2. Sandwich with FGM Core

In this example, we analyze the bending of a square sandwich B-type plate with thicknessh. The bottom skin is aluminum (
E
m
=
70
G
P
a
) with thickness h
b
=
0.1
h and the top skin is alumina (
E
c
=
380
G
P
a
) with thickness h
t
=
0.1
h. The core is in FGM with volume fraction of the ceramic according to (2).The transverse displacement and the normal stresses are presented in normalized form as(36)
τ
¯
x
y
=
-
h
2
q
a
2
τ
x
y
x
=
a
,
y
=
a
,
z
=
0
,
z
=
h
/
2
,
τ
¯
z
x
=
h
q
a
τ
z
x
x
=
a
,
y
=
a
,
z
=
-
h
/
N
e
/
2
.The convergence study was performed forσ
x
z and transverse displacement w
(
0
) considering p
=
4, a
/
h
=
100, and the order K of the plate theory of 1, 3, 5, and 7. Results are presented in Tables 4 and 5. We consider that the 5-order plate theory should be used in the following computation.Table 4
w convergence study for the bending analysis of B-type plate using higher-order plate theory, p
=
4, and a
/
h
=
100.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (K
=
5
/
6) [37]
Ref. [39]
Ref. [41]
w
0.7796
0.7797
0.7797
0.7797
0.7796
0.7797
0.7797Table 5
σ
x
z convergence study for the bending analysis of B-type plate using higher-order plate theory, p
=
4, and a
/
h
=
100.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (K
=
5
/
6) [37]
Ref. [39]
Ref. [41]
σ
x
z
0.1176
0.1949
0.1862
0.1819
0.1877
0.2398
0.2432In Table6 we present the values of σ
x
z and out-of-plane displacement for various values of exponent p of the material power-law (p
=
0,0.5,1
,
4,10) and various thickness to side ratios (a
/
h
=
4,10,100) according to the present higher-order theory using the 5-order plate theory. Results are tabulated and compared with available references. It demonstrates the results are in excellent agreement with the meshfree results [44] and other reference values under various thickness to span ratio (a
/
h) and four and five different kinds of parameters (p).Table 6
Square B-type plate in bending with consideringσ
z
x and w use the present higher-order plate theory.
p
a
/
h
σ
z
x
(
h
/
6
)
w
(
0
)
4
10
100
4
10
100
0
Meshfree [37]
0.2208
0.2227
0.2228
0.4447
0.3711
0.3568
Present
0.2392
0.2411
0.2415
0.4655
0.3748
0.3579
0.5
Meshfree [37]
0.2546
0.2581
0.2585
0.6168
0.5238
0.5058
Present
0.2406
0.2417
0.2418
0.6409
0.5240
0.5023
1
Ref. [41]
0.2613
0.2605
0.2603
0.7628
0.6324
0.6072
CLPT [37]
0.0000
0.0000
0.0000
0.6070
0.6070
0.6070
FSDT (k
=
5
/
6) [37]
0.2458
0.2458
0.2458
0.7738
0.6337
0.6073
Ref. [39] N
=
4
0.2604
0.2594
0.2593
0.7628
0.6324
0.6072
Ref. [40]
0.2742
0.2788
0.2793
0.7416
0.6305
0.6092
Meshfree [37]
0.2745
0.2789
0.2795
0.7417
0.6305
0.6092
Present
0.2287
0.2290
0.2290
0.7766
0.6339
0.6073
4
Ref. [41]
0.2429
0.2431
0.2432
1.0934
0.8321
0.7797
CLPT [37]
0.0000
0.0000
0.0000
0.7792
0.7792
0.7792
FSDT (k
=
5
/
6) [37]
0.1877
0.1877
0.1877
1.0285
0.8191
0.7796
Ref. [39] N
=
4
0.2400
0.2398
0.2398
1.0930
0.8307
0.7797
Ref. [40]
0.2723
0.2778
0.2785
1.0391
0.8202
0.7784
Meshfree [37]
0.2696
0.2747
0.2753
1.0371
0.8199
0.7784
Present
0.1850
0.1861
0.1862
1.1034
0.8311
0.7797
10
Ref. [41]
0.2150
0.2174
0.2179
1.2232
0.8753
0.8077
CLPT [37]
0.0000
0.0000
0.0000
0.8070
0.8070
0.8070
FSDT (k
=
5
/
6) [37]
0.1234
0.1234
0.1234
1.1109
0.8556
0.8075
Ref. [39] N
=
4
0.1932
0.1944
0.1946
1.2172
0.8740
0.8077
Ref. [40]
0.2016
0.2059
0.2064
1.1780
0.8650
0.8050
Meshfree [37]
0.1995
0.2034
0.2039
1.1752
0.8645
0.8050
Present
0.1835
0.1862
0.1866
1.2372
0.8755
0.8077In Figure5, we also present the evolution of the displacement and stresses across the thickness direction according to present shear deformation theory for various values of the exponent p of a plate with side-to-thickness ratio a
/
h
=
100, using the 5-order plate theory. Compared with the meshfree method, all the values are consisted with each other except the late one. The major difference is located in internal boundary (h
=
2
h
/
5) on out-of-plane stress σ
z
x. In addition, the present solution of the stresses σ
z
x is very close to zero at upper skin and bottom skin (p
=
1,4
,
10), and the meshfree solution is away from the origin of axis. In conclusion, these results indicate that the OHSNDPT is a good choice to deal with the FGM plate problem.Figure 5
Square B-type plate subjected to sinusoidal load at the top, witha
/
h
=
100. Dimensionless displacement (w
¯) through the thickness direction according to present higher-order theory for different values of p.


#### 5.1.3. Isotropic FGM Plate Compared with Laminated Plate

To demonstrate the availability of the present orthogonal higher-order shear and deformation plate theory, more cases of the modulus form are considered, linear model, exponential model, and reciprocal model (refer to [47]).The attention is focused on the influence of the three distributions of the elastic modulus on the displacement and stress fields in the plate. Young’s moduli at the upper and lower surfaces of the plate are given asE
(
0
)
=
1
G
P
a and E
(
h
)
=
10
G
P
a. Accordingly, material parameters for FGM models used in the present study are determined as follows:(a)
Liner modelE
(
z
)
=
E
0
(
p
(
z
/
h
+
1
/
2
)
+
1
), E
0
=
E
(
0
)
=
1
G
P
a, p
=
9
(b)
Exponential model,E
(
z
)
=
E
0
e
p
∗
(
z
/
h
+
1
/
2
), E
0
=
E
(
0
)
=
1
G
P
a, p
=
ln
⁡
E
h
/
E
0
=
ln
⁡
10
(c)
Reciprocal model,E
(
z
)
=
E
0
/
p
∗
(
z
/
h
+
1
/
2
)
+
1, E
0
=
E
(
0
)
=
1
G
P
a, p
=
E
(
0
)
/
E
(
h
)
-
1
=
-
0.9Based on the above three models of Young’s modulus, the variation of displacements u and w, stresses σ
x
,
τ
x
y
,
σ
z, and τ
z
x at a chosen location (
x
/
a
=
1
/
4
,
y
/
b
=
1
/
4
), as a function of the z-coordinate, are shown in Figure 6 for the plate (
h
/
a
=
0.1
).Figure 6
D-type square plate for several exponents for the powerp compared with the type E schemes.
(a)
(b)
(c)
(d)From Figure6, the following observations can be made:(I)
The magnitude of the vertical displacement (deflection) for the reciprocal model is biggest while that of the linear model is the least (Figure6(a)). This reveals that the bending rigidities of the plate are different for these three graded models. The linear model is most rigid and the reciprocal model is softest.
(II)
The in-plane stresses,σ
x
x (Figure 6(c)) and τ
x
y (Figure 6(d)), are essentially nonlinear over the thickness of the plate for all these three models of functionally graded materials. The in-plane stress concentrations in the plate are quite different for different grade models of Young’s modulus. For example, the maximum compressive stress for the reciprocal model is about 1.8 times that of the linear model (see Figure 6(c), especially 63.0921 and 34.2402 at the center point). This enables an optimal design of the plate by selecting appropriately graded models of materials.
(III)
Compared with the 3D closed-form solutions [47], the present solutions of the OHOSNDPT are well coincident with the reference data from Figure 6. Not only the in-plane stresses, but also σ
x and τ
x
y reach an identical with the 3D closed-form solutions, but also the out-of-plane stress, τ
z
x, is in very good agreement with the similar values in the literature.In short, the present approach can actually reflect the distribution situation of displacements and stresses. In addition, it is also in very good agreement with similar theories in the literature. Thus, the present method is a good idea worth trying.


#### 5.1.4. Power-Law Fitted with the Straight Line

In this section, we state another scheme to illustrate the superiority of the OHOSNDPT.Because any function can be fitted with straight lines, the stiffness functionψ
(
z
) of the FGM plate with arbitrary form equivalents to a piecewise function ψ
¯
(
z
) with the function is a linear function ψ
¯
k
(
z
) in each interval. Moreover, the present orthogonal higher-order shear and nonamial plate theory also be suitable for the piecewise function ψ
¯
(
z
) in laminated FGM plates.As shown in Figure7, the variations of Young’s modulus are along the thickness of the power-law FGM plates although their values at the upper and lower surfaces of the plate are fixed and the corresponding piecewise linear functions ψ
¯
(
z
) are chosen.Figure 7
D-type square plate for several exponents for the powerp compared with the type E schemes.Constant coefficients are usually used in the each layer when the FGM plate can be seen as the laminated plates [47]. This scheme facilitates settling the FGM models, but some shortcomings are also exposed in Figure 8. Figure 8 shows the variation of the longitudinal displacements w, in-plane stresses σ
x
x and σ
x
y and out-plane stresses σ
y
z of the Tape A simple supported FGM square plate (k
=
3, p
=
10 and a
/
h
=
10).(i)
The constant coefficients of the FGM material properties in each layer are generally taken to be the value of the midpoint of the layer interval. For the class of material which the volume fraction varies steeply (p
=
10), the errors between the laminated plate models and the FGM properties obviously exist. Thus, the errors between the scheme of straight line and the analytical solution are smaller than the constant coefficients ones and it is verified from Figure 8(a).
(ii)
If the constant coefficients of the laminated plates model are used, the paramount material properties, the expression of Young’s modulus, and the density function are not continuous functions on vertical directionz. The displacement model is usually supposed as the continuous function on z, the discontinuous of stresses are obtained according to geometric equations and physics equations. Obviously, this is contrary to reality. From Figures 8(b) and 8(c), it can be observed that the scheme of straight line can match the analytical solution. The curves of the constant coefficients appear as interrupted phenomena on h
/
6, and the error enlarges with the increasing of z. On transverse stress σ
y
z (Figure 8(d)), the analytical solution and the schemes of the straight lines and the constant coefficients all occur the discontinuous situation on h
/
6, but the performance of the straight lines is better than the constant coefficients and close to the analytical solutions.Figure 8
D-type square plate for several exponents for the powerp compared with the type E schemes.
(a)
(b)
(c)
(d)From Table7, some phenomena can be obtained:(i)
Whenk is small (k
=
2,3
,
5), the scheme of the straight lines is better than the approach of the constant coefficients.
(ii)
Whenk is large (k
=
200,400), all two schemes achieve the convergence solutions.
(iii)
The two scheme coupled with OHOSNDPT are superior to the common scheme of the constant coefficients, while the familiar approach takesk
=
1000 to get the stable values.Table 7
The convergence study for the laywise number of the B-type plate using higher-order plate theory,p
=
10, and a
/
h
=
4.
Number
2
3
5
10
50
100
200
400
Analytical
Laywise (Error)
1.8480
1.6586
1.5094
1.4311
1.4038
1.4029
1.4027
1.4026
1.4026
31.76%
18.25%
7.61%
2.03%
0.09%
0.02%
0.01%
0.00%
Linear (Error)
1.1328
1.2422
1.3206
1.3770
1.4015
1.4023
1.4025
1.4026
19.24%
11.44%
5.85%
1.83%
0.08%
0.02%
0.01%
0.00%We draw a conclusion from the above example, the OHOSNDPT is a very suitable FGM plate theory.


### 5.2. Free Vibrations

In this example, we study the free vibration of a simply supported isotropic FGM square plate (a
=
b
=
1) of type A. The plate is graded from aluminum (bottom) to zirconia/alumina (top). The values for metal and ceramics in the FGM plate are listed in Table 8.Table 8
Material properties of the used FGM plate.
Material
Properties
E
(
G
P
a
)
ν
ρ
(
k
g
/
m
3
)
Aluminum (A
l)
70
0.3
2702
Alumina (A
l
2
O
3)
380
0.3
3800
Zirconia (Z
r
O
2)
200
0.3
5700We consider the Mori–Tanaka homogenization scheme (36), as in Vel and Batra [51] (here considered to be the exact solution), and as in Neves et al. [40] and Qian et al. [44].The frequencyw has been nondimensionalized as follows:(37)
w
¯
=
w
h
ρ
m
E
m
.


#### 5.2.1. Isotropic FGM Plate

In Table9 we present the results obtained with the theories considered and different values of p for a side-to-thickness ratio a
/
h
=
5. Results show that the convergence values are obtained with K
=
3. It is noteworthy that the present solutions are almost the same as values of [48] for p
=
1 regardless of the ratio length-thickness.Table 9
Fundamental frequency of a square simply supported A-type square plate (Al/ZrO2) using the present higher-order theory.
Mode
p
=
1
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
p
=
2
p
=
3
p
=
5
Exact 3D [48]
0.0153
0.0596
0.2192
0.2197
0.2211
0.2225
HOSNT-15 [49]
0.0154
0.0596
0.2191
0.2196
0.2211
0.2225
Ref. [37]
0.0153
0.0596
0.2193
0.2200
0.2215
0.2230
Ref. [44]
0.0149
0.0584
0.2152
0.2153
0.2172
0.2194
Ref. [42]
0.0158
0.0619
0.2276
0.2264
0.2276
0.2291
Present (K
=
1)
0.0159
0.0621
0.2307
0.2295
0.2309
0.2326
Present (K
=
3)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2272
Present (K
=
5)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2271
Present (K
=
7)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2271To further compare the exact HSDT solution and the present OHSNDPT, we introduce more reliable data. In Table10, based on the present closed-form solutions, exact HSDT solution with refined plate theory and finite element method, numerical results have been performed for simple supported Al/Al2O3 square plates when p
=
1. A well-known commercially available FEM package was used to investigate 3D free vibration of FGM square plates in [48].Table 10
The natural frequency of A-type square simply supported plate (A
l
/
A
l
2
O
3) using the higher-order theory.
a
/
h
Model
1
(
1,1
)
2
(
1,2
)
3
(
2,1
)
4
(
2,2
)
5
(
1,3
)
6
(
3,1
)
7
(
2,3
)
8
(
3,2
)
9
(
1,4
)
1
0
(
4,1
)
20
Present
359.90
889.22
889.22
1406.6
1745.3
1745.3
2244.4
2244.4
2894.3
2894.3
Ref. [43]
359.97
889.39
889.39
1406.9
1745.7
1745.7
2244.9
2244.9
—
—
FEM [43]
357.37
883.58
883.58
1398.6
1736.2
1736.2
2234.1
2234.1
—
—
Error
-
0.71
%
-
0.64
%
-
0.64
%
-
0.57
%
-
0.52
%
-
0.52
%
-
0.46
%
-
0.46
%
10
Present
703.30
1685.4
1685.4
2596.2
3170.5
3170.5
3990.0
3990.0
5016.1
5016.1
Ref. [43]
703.45
1685.7
1685.7
2596.6
3170.7
3170.7
3989.9
—
—
—
FEM [43]
699.30
1679.7
1679.7
2592.5
3169.6
3169.6
3994.8
—
—
—
Error
-
0.57
%
-
0.34
%
-
0.34
%
-
0.14
%
-
0
.
0
3
%
-
0.03
%
0.12%
5
Present
1298.1
2871.7
2871.7
4186.8
4969.4
4969.4
6041.0
6041.0
7327.4
7327.4
Ref. [43]
1298.3
2870.5
2870.5
4181.8
4960.5
4960.5
—
—
—
—
FEM [43]
1296.3
2883.0
2883.0
4216.1
5010.6
5010.6
—
—
—
—
Error
-
0.14
%
0.39%
0.39%
0.69%
0.82%
0.82%The length of square plates is 1 m. Three different thicknesses 0.05 m (corresponding to thin plates), 0.1 m and 0.2 m (corresponding to moderately thick plates) have been used. The calculations of the present scheme are obtained for first ten natural frequencies. “—” values in Table10 indicate no calculation in corresponding reference. The percentage difference given in Table 10 is defined as follows: E
r
r
o
r
(
%
)
=
[
F
E
M
-
(
s
o
l
u
t
i
o
n
)
]
/
F
E
M.An excellent agreement is observed between the present plate theory, exact HSDT solution and the FEM. It is seen that the exact results of proposed method are close to the FEM analysis results at lower and higher frequencies. Frequencies rise with an increase in the thickness of plate. This phenomena originates from the increasing the rigidity of plate. Frequencies decrease when less restraining boundary is used at the edges of square plates. This is due to the fact that higher constraints at the edges increase the flexural rigidity of the plate, leading to a higher frequency response. Compared with the exact HSDT solution, the present closed-form solution is closer to the FEM solution.


#### 5.2.2. Isotropic FGM Plate Compared with Laminated Plate

In this subsection the higher-order plate theories are used for the free vibration analysis of simply supported functionally of three models of type D with side-to-thickness ratioa
/
h
=
10. All the material properties is the same as that in Section 5.1.3 and the density function is supposed with ρ
(
z
)
=
ρ
0
ψ
(
z
), that is, ρ
(
h
/
2
)
=
10
∗
ρ
(
-
h
/
2
).In Table11, nondimensional first 10 frequencies of plate are computed for three different parameter a
/
h
=
20,10,5 in the Tape D. Some views can be noticeable:(i)
The order of the first ten frequency of the three model all comply with the order of isotropic square plate (m
2
+
n
2). It shows the variety of the properties of material on z has nothing to do with the order of the frequencies.
(ii)
Nondimensional results are hardly affected bya
/
h.
(iii)
Unlike static problem, the longitudinal displacementw for the reciprocal model is about 1.8 times that of the linear model. The frequencies of them are very similar to each other for fixed (
m
,
n
). It is because the ratio of the upper material and the bottom material of the density function is set as 10, which is consistent with the ratio of elastic modulus.Table 11
First ten frequencies of type D square plate with variousa
/
h using the higher-order theory.
Source
1
(
1,1
)
2
(
1,2
)
3
(
2,1
)
4
(
2,2
)
5
(
1,3
)
6
(
3,1
)
7
(
2,3
)
8
(
3,2
)
9
(
1,4
)
1
0
(
4,1
)
Linear model
a
/
h
=
20
0.8481
2.0986
2.0986
3.3243
4.1284
4.1284
5.3158
5.3158
6.8659
6.8659
a
/
h
=
10
0.8311
2.0012
2.0012
3.0947
3.7874
3.7874
4.7796
4.7796
6.0269
6.0269
a
/
h
=
5
0.7737
1.7284
1.7284
2.5333
3.0136
3.0136
3.6719
3.6719
4.4619
4.4619
Exponent model
a
/
h
=
20
0.8504
2.1025
2.1025
3.3280
4.1310
4.1310
5.3155
5.3155
6.8596
6.8596
a
/
h
=
10
0.8320
1.9982
1.9982
3.0833
3.7688
3.7688
4.7485
4.7485
5.9769
5.9769
a
/
h
=
5
0.7708
1.7121
1.7121
2.5011
2.9708
2.9708
3.6139
3.6139
4.3853
4.3853
Reciprocal model
a
/
h
=
20
0.9214
2.2722
2.2722
3.5877
4.4464
4.4464
5.7559
5.7559
7.3460
7.3460
a
/
h
=
10
0.8969
2.1357
2.1357
3.2725
3.9843
3.9843
4.9939
4.9939
6.2496
6.2496
a
/
h
=
5
0.8181
1.7836
1.7836
2.5767
3.0444
3.0444
3.6811
3.6811
4.4410
4.4410


### 5.3. Buckling

#### 5.3.1. Sandwich Plate with FGM Skins

In the next examples the higher-order plate theories are used for the buckling analysis of simply supported functionally graded sandwich square plates (a
=
b
=
1) of type C with side-to-thickness ratio a
/
h
=
10. The uni-and biaxial critical buckling loads are analsed.The material properties areE
m
=
70
E
0 (aluminum) for the metal and E
c
=
380
E
0 (alumina) for the ceramic being E
0
=
1
G
P
a. The law-of-mixtures ((19)–(21)) was used for the computation of Young’s modulus for each layer. The nondimensional parameter used is(38)
P
¯
=
P
a
2
100
h
3
E
0
.An initial convergence study with the higher-order theory was conducted for each buckling load type considering the orderK of the plate theory is 1, 3, 5, and 7. The uniaxial case is presented in Table 12 for the 2
-
2
-
1 sandwich with p
=
5 and the biaxial case is presented in Table 13 for the 1
-
2
-
1 sandwich with p
=
1. Further results are obtained by considering K equal to 5, which seems acceptable by the convergence study.Table 12
Convergence study for the uniaxial buckling load of a simply supported 2-2-1 sandwich square plate with FGM skins andp
=
5 case using the higher-order theory.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT [37]
TSDPT [45]
SSDPT [46]
P
-
4.11230
4.10740
4.10605
4.10449
4.09285
4.11209
4.11280Table 13
Convergence study for the uniaxial buckling load of a simply supported 1-2-1 sandwich square plate with FGM skins andp
=
1 case using the higher-order theory.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT [37]
TSDPT [45]
SSDPT [46]
P
-
3.76605
3.75329
3.75304
3.75290
3.74182
3.75328
3.75314The critical buckling loads obtained from the present approach are tabulated in Tables14 and 15 for various power-law exponents p and thickness ratios. Both tables include results obtained from classical plate theory (CLPT), first-order shear deformation plate theory (FSDPT, the shear correction factor is taken as 5
/
6), Reddy’s higher-order shear deformation plate theory (TSDPT) [10], and Zenkour’s sinusoidal shear deformation plate theory (SSDPT) [29]. Table 14 refers to the uniaxial buckling load and Table 15 refers to the biaxial buckling load.Table 14
Uniaxial buckling load of simply supported plate of C-type using the 5-order theory.
p
Theory
P
-
1-0-1
2-1-2
2-1-1
1-1-1
2-2-1
1-2-1
0
CLPT [37]
13.73791
13.73791
13.73791
13.73791
13.73791
13.73791
FSDPT [37]
13.00449
13.00449
13.00449
13.00449
13.00449
13.00449
TSDPT [45]
13.00495
13.00495
13.00495
13.00495
13.00495
13.00495
SSDPT [46]
13.00606
13.00606
13.00606
13.00606
13.00606
13.00606
Meshfree [37]
13.95287
13.95287
13.95287
13.95287
13.95287
13.95287
13.00516
13.00516
13.00516
13.00516
13.00516
13.00516
0.5
CLPT [37]
7.65398
8.25597
8.56223
8.78063
9.18254
9.61525
FSDPT [37]
7.33732
7.91320
8.20015
8.41034
8.78673
9.19517
TSDPT [45]
7.36437
7.94084
8.22470
8.43645
8.80997
9.21681
SSDPT [46]
7.36568
7.94195
8.22538
8.43712
8.81037
9.21670
Meshfree [37]
7.16207
7.71627
7.98956
8.22133
8.55172
8.94190
7.36361
7.94047
8.22048
8.43637
8.80980
9.21651
1
CLPT [37]
5.33248
6.02733
6.40391
6.68150
7.19663
7.78406
FSDPT [37]
5.14236
5.81379
6.17020
6.43892
6.92571
7.48365
TSDPT [45]
5.16713
5.84006
6.19394
6.46474
6.94944
7.50656
SSDPT [46]
5.16846
5.84119
6.19461
6.46539
6.94980
7.50629
Meshfree [37]
5.06137
5.71135
6.05467
6.31500
6.78405
7.31995
5.16489
5.83877
6.19209
6.48150
6.94876
7.50579
5
CLPT [37]
2.73080
3.10704
3.48418
3.65732
4.21238
4.85717
FSDPT [37]
2.63842
3.02252
3.38538
3.55958
4.09285
4.71475
TSDPT [45]
2.65821
3.04257
3.40351
3.57956
4.11209
4.73469
SSDPT [46]
2.66006
3.04406
3.40449
3.58063
4.11288
4.73488
Meshfree [37]
2.63652
3.00791
3.36225
3.53005
4.05070
4.64701
2.64204
3.02902
3.38599
3.57260
4.10499
4.73116
10
CLPT [37]
2.56985
2.80340
3.16427
3.25924
3.79238
4.38221
FSDPT [37]
2.46904
2.72626
3.07428
3.17521
3.68890
4.26040
TSDPT [45]
2.48727
2.74632
3.09190
3.19471
3.70752
4.27991
SSDPT [46]
2.48928
2.74844
3.13443
3.19456
3.14574
4.38175
Meshfree [37]
2.47216
2.72046
3.06067
3.15761
3.66166
4.20550
2.47035
2.72423
3.06339
3.18364
3.09531
4.27493Table 15
Biaxial buckling load of simply supported plate of C-type using the 5-order theory.
p
Theory
P
-
Source
1-0-1
2-1-2
2-1-1
1-1-1
2-2-1
1-2-1
0
CLPT [37]
6.86896
6.86896
6.86896
6.86896
6.86896
6.86896
FSDPT [37]
6.50224
6.50224
6.50224
6.50224
6.50224
6.50224
TSDPT [45]
6.50248
6.50248
6.50248
6.50248
6.50248
6.50248
SSDPT [46]
6.50303
6.50303
6.50303
6.50303
6.50303
6.50303
Meshfree [37]
6.47643
6.47643
6.47643
6.47643
6.47643
6.47643
6.50258
6.50258
6.50258
6.50258
6.50258
6.50258
0.5
CLPT [37]
3.82699
4.12798
4.28112
4.39032
4.59127
4.80762
FSDPT [37]
3.66866
3.95660
4.10007
4.20517
4.39336
4.59758
TSDPT [45]
3.68219
3.97042
4.11235
4.21823
4.40499
4.60841
SSDPT [46]
3.68284
3.97097
4.11269
4.21856
4.40519
4.60835
Meshfree [37]
3.58104
3.85813
3.99478
4.09639
4.27586
4.47095
3.68180
3.97024
4.11204
4.21819
4.40490
4.60826
1
CLPT [37]
2.66624
3.01366
3.20195
3.34075
3.59831
3.89203
FSDPT [37]
2.57118
2.90690
3.08510
3.21946
3.46286
3.74182
TSDPT [45]
2.58357
2.92003
3.09697
3.23237
3.47472
3.75328
SSDPT [46]
2.58423
2.92060
3.09731
3.23270
3.47490
3.75314
Meshfree [37]
2.53069
2.85568
3.02733
3.15750
3.39202
3.65998
2.58245
2.91939
3.09605
3.23212
3.47438
3.75290
5
CLPT [37]
1.36540
1.55352
1.74209
1.82866
2.10619
2.42859
FSDPT [37]
1.31921
1.51126
1.69269
1.77979
2.04642
2.35737
TSDPT [45]
1.32910
1.52129
1.70176
1.78978
2.05605
2.36734
SSDPT [46]
1.32910
1.52203
1.70224
1.79032
2.05644
2.36744
Meshfree [37]
1.31826
1.50395
1.68128
1.76502
2.02535
2.32351
1.32102
1.51451
1.69299
1.78630
2.05225
2.36558
10
CLPT [37]
1.28493
1.40170
1.58214
1.62962
1.89619
2.19111
FSDPT [37]
1.23452
1.36313
1.53714
1.58760
1.84445
2.13020
TSDPT [45]
1.24363
1.37316
1.54595
1.59736
1.85376
2.13995
SSDPT [46]
1.24475
1.37422
1.56721
1.59728
1.57287
2.19087
Meshfree [37]
1.23608
1.36023
1.53034
1.57880
1.83083
2.10275
1.23518
1.36211
1.53170
1.59182
1.84765
2.13746A good agreement between the present solution and references considered; especially [10, 29] is obtained. This allows us to conclude that the present higher-order plate theory is good for the modeling of simply supported sandwich FGM plates with OHOSNDPT as a good formulation. This study also lead us to conclude that the thickness stretching effect has a strong influence on the buckling analysis of sandwich FGM plates as gives higher fundamental buckling loads.The isotropic fully ceramic plate (first line on Tables14 and 15) has the higher fundamental buckling loads. As the core thickness to the total thickness of the plate ratio (
(
h
2
-
h
1
)
/
h
) increases, the buckling loads increase as well. Considering each column of both tables we may conclude that the critical buckling loads decrease as the power-law exponent p increases. By comparing Tables 14 and 15 we also conclude that the biaxial buckling load of simply supported sandwich square plate with FGM skins is half the uniaxial one for the same plate.


#### 5.3.2. Isotropic FGM Plate Compared with Laminated Plate

In this subsection, the higher-order plate theories are used for the buckling analysis of simply supported functionally of two models of type D and E with side-to-thickness ratioa
/
h
=
10. All the material properties are taken as in Section 5.1.3.The convergence study with the higher-order theory was also conducted for each buckling load type considering the orderK of the plate theory from 1 to 7. The uniaxial case is presented in Table 16 for the Linear model with a
/
h
=
20,10,5. Further results are obtained by considering K equals 3, which seems acceptable by the convergence study.Table 16
Convergence study for the uniaxial buckling load of a simply supported square plate with linear model using the higher-order theory.
Scheme
K
=
1
K
=
2
K
=
3
K
=
4
K
=
5
K
=
6
K
=
7
a
/
h
=
20
30.6144
30.5958
30.5585
30.5584
30.5584
30.5584
30.5584
a
/
h
=
10
14.9028
14.8679
14.7979
14.7978
14.7978
14.7977
14.7977
a
/
h
=
5
6.7392
6.6849
6.5739
6.5739
6.5738
6.5738
6.5738In Table17, nondimensional uniaxial buckling load and biaxial buckling load of plate are computed for various a
/
h of three different models.Table 17
Bi and iniaxial buckling load of simply supported plate of type D using the higher-order theory.
Model
Reciprocal
Exponent
Linear
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
Biaxial
15.2792
7.3989
3.2869
10.9162
5.2696
2.3170
8.3944
4.0156
1.7136
Uniaxial
30.5585
14.7979
6.5739
21.8325
10.5393
4.6339
16.7888
8.0313
3.4273It is well known that as the plate thickness increases, the shear effect which leads to reduction of the critical buckling stress increases. It can be seen from Tables16 and 17 that although increasing the plate thickness decreases P
¯
x in all models, it can be concluded that for compressive biaxial and uniaxial loadings consideration of the prebuckling deformations in the stability equations increases the shear effect.Another interesting phenomenon appears between the static problem and the buckling situation. Whena
/
h
=
10, that is, the FGM is a thin plate, the rate of the buckling stresses between the reciprocal model and the linear model (the value is 1.8425) is almost the same as the proportion of the displacement at the center point (the more accurate value 1.8426 is than the number in Section 5.1.4) at the same K. When a
/
h
=
5, that is, the FGM is a moderately thick plate, the former ratio is 1.9181 and the latter one becomes 1.9198 (4.6322/2.4128). In summary, it not only illustrate that the present OHOSNDPT is high accuracy both the static problem and the buckling situation, but also reflect the stability for the thin plate and moderately thick plate.


## 5.1. Bending

In the following static examples, we consider that the plate is subjected to a bisinusoidal transverse mechanical load of amplitude loadq
t
x
,
y
=
q
¯
t
sin
⁡
π
x
/
a
sin
⁡
π
y
/
b applied at the top of the plate with q
¯
t
=
1. It should be noted that the load is applied at the top surface (
z
=
h
/
2
) as shown in Figure 2.Due to the symmetry of the problem, the displacementv, stresses σ
y, and transverse shear stresses τ
z
y are not described since their distributions along the plate thickness direction are similar to those of u
,
σ
x and τ
z
x, respectively. Considering σ
x is the first-order derivative function of the displacement u, only the stress σ
x should be discussed in this paper.


### 5.1.1. Isotropic FGM Square Plate

In this example, an isotropic FGM square plate of type-A is considered. The plate is graded from aluminumE
m
=
70
G
P
a at the bottom to alumina E
c
=
380
G
P
a at the top. The transverse displacement, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as(35)
w
¯
=
10
h
3
E
c
q
a
4
w
,
σ
¯
x
x
=
h
2
q
a
2
σ
x
x
.An initial convergence study was performed forσ
x
x
h
/
3 and transverse displacement w
(
0
) at the center of the plate, considering p
=
1, a
/
h
=
10, and the order K of the plate theory of 1, 3, 5, and 7. Compared with the references [37] and [39], results are presented in Tables 1 and 2. As seen in these tables, it is sufficient to use K
=
5 (we consider that K
=
5 should be used in the following computation in this section).Table 1
w
- convergence study for the bending analysis of A-type plate using higher-order plate theory, p
=
1, and a
/
h
=
10.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT(
k
=
5
/
6
) [37]
[39]
Meshfree [37]
w
-
0.5846
0.5891
0.5891
0.5891
0.5889
0.5890
0.5868Table 2
σ
-
x
x convergence study for the bending analysis of A-type plate using higher-order plate theory, p
=
1, and a
/
h
=
10.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (k
=
5
/
6) [37]
[39]
Meshfree [37]
σ
-
x
x
1.4968
1.4904
1.4908
1.4908
2.0150
1.5064
1.4917In Table3, we present results for σ
¯
x
x and transverse displacement for various exponents p of the power-law (18) considering the order of the plate theory equal 7. The considered side-to-thickness ratios (
a
/
h
) are 4, 10, and 100, which means thickness h equals 0.25, 0.1, and 0.01, respectively. Results are compared with the classical plate theory (CLPT), the first-order shear deformation theory (FSDT) with a correction factor k
=
5
/
6, and those from Zenkour’s generalized shear deformation theory [38], and those from Carrera et al. [25, 39], and Neves et al. [40].Table 3
A-type plate in bending. Theσ
z
x and w according to present higher-order plate theory.
p
a
/
h
σ
x
x
(
h
/
3
)
w
(
0
)
4
10
100
4
10
100
0
Meshfree [37]
≠0
0.5278
1.3176
13.161
0.3665
0.2942
0.2803
Present
0.5144
1.3123
13.1713
0.3800
0.2961
0.2804
0.5
Meshfree [37]
≠0
0.5860
1.4680
14.673
0.5493
0.4548
0.4365
Present
≠0
0.5724
1.4591
14.6446
0.5678
0.4538
0.4326
1
Ref. [41]
≠0
0.6221
14.969
14.969
0.7171
0.5875
0.5625
CLPT [37]
0
0.8060
2.0150
20.150
0.5623
0.5623
0.5623
FSDT (k
=
5
/
6) [37]
0
0.8060
2.0150
20.150
0.7291
0.5889
0.5625
Ref. [39] N
=
4
≠0
0.6221
1.5064
14.969
0.7171
0.5875
0.5625
Ref. [40]
≠0
0.5925
1.4945
14.969
0.6997
0.5845
0.5624
Meshfree [37]
≠0
0.5911
1.4917
14.945
0.7020
0.5868
0.5647
Present
0.5841
1.4908
14.9677
0.7318
0.5891
0.5625
4
Ref. [41]
≠0
0.4877
1.1971
11.923
1.1585
0.8821
0.8286
CLPT [37]
0
0.6420
1.6049
16.049
0.8281
0.8281
0.8281
FSDT (k
=
5
/
6) [37]
0
0.6420
1.6049
16.049
1.1125
0.8736
0.828
Ref. [39] N
=
4
≠0
0.4877
1.1971
11.923
1.1585
0.8821
0.8286
Ref. [40]
≠0
0.2723
0.2778
0.2785
1.0391
0.8202
0.7784
Meshfree [37]
≠0
0.4330
1.1588
11.737
1.1108
0.8700
0.8240
Present
0.4451
1.1787
11.9209
1.1732
0.8831
0.8287
10
Ref. [41]
≠0
0.3695
0.8965
8.9077
1.3745
1.0072
0.9361
CLPT [37]
0
0.4796
1.1990
11.990
0.9354
0.9354
0.9354
FSDT (k
=
5
/
6) [37]
0
0.4796
1.1990
11.990
1.3178
0.9966
0.9360
Ref. [39] N
=
4
≠0
0.1478
0.8965
8.9077
1.3745
1.0072
0.9361
Ref. [40]
≠0
0.3227
1.1783
11.932
1.3490
0.8750
0.8286
Meshfree [37]
≠0
0.3097
0.8462
8.6010
1.3334
0.9888
0.9227
Present
≠0
0.3238
0.8766
8.9058
1.4026
1.0094
0.9362Date shown in Table3 clearly illustrated, a subtle difference appears in the moderately thick plate between them and the present higher-order plate theory are in good agreement with those from the Refs. When the side-to-thickness ratio is 100, the values, whether it is on transverse displacements or stresses, are almost the same as the values of [50].In Figure4, we present the evolution of the displacement and stresses across the thickness direction according to present shear deformation theory for various values of the exponent p and side-to-thickness ratio a
/
h
=
4, using the 5-order plate theory. Compared with the results of the meshfree method in [44], it shows that the values of the in-plane stress σ
x
x, σ
x
y and out-plane stress σ
z
x are consistent with the reference value. Strangeness is the presence of a significant difference between the present theory (K
=
5) and the meshless solutions on the longitudinal displacement w on p
=
10. From Figure 4, we find out the relative error of displacement w can up to 9% and the maximum relative error of stresses between them is less than 3%. This discord phenomenon is due to the computational error and the instability of the numerical method. On the other hand, it demonstrates the analytical solution of the orthogonal higher-order shear deformable plate theory is necessary and feasible.Figure 4
A-type square plate subjected to sinusoidal load at the top, witha
/
h
=
4, nondimensional displacement (w
¯) and stresses through the thickness to present OHOSNDPT for different values of p.In summary, it can be concluded that the present higher-order theory provides excellent solution for the isotropic FGM plates.


### 5.1.2. Sandwich with FGM Core

In this example, we analyze the bending of a square sandwich B-type plate with thicknessh. The bottom skin is aluminum (
E
m
=
70
G
P
a
) with thickness h
b
=
0.1
h and the top skin is alumina (
E
c
=
380
G
P
a
) with thickness h
t
=
0.1
h. The core is in FGM with volume fraction of the ceramic according to (2).The transverse displacement and the normal stresses are presented in normalized form as(36)
τ
¯
x
y
=
-
h
2
q
a
2
τ
x
y
x
=
a
,
y
=
a
,
z
=
0
,
z
=
h
/
2
,
τ
¯
z
x
=
h
q
a
τ
z
x
x
=
a
,
y
=
a
,
z
=
-
h
/
N
e
/
2
.The convergence study was performed forσ
x
z and transverse displacement w
(
0
) considering p
=
4, a
/
h
=
100, and the order K of the plate theory of 1, 3, 5, and 7. Results are presented in Tables 4 and 5. We consider that the 5-order plate theory should be used in the following computation.Table 4
w convergence study for the bending analysis of B-type plate using higher-order plate theory, p
=
4, and a
/
h
=
100.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (K
=
5
/
6) [37]
Ref. [39]
Ref. [41]
w
0.7796
0.7797
0.7797
0.7797
0.7796
0.7797
0.7797Table 5
σ
x
z convergence study for the bending analysis of B-type plate using higher-order plate theory, p
=
4, and a
/
h
=
100.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (K
=
5
/
6) [37]
Ref. [39]
Ref. [41]
σ
x
z
0.1176
0.1949
0.1862
0.1819
0.1877
0.2398
0.2432In Table6 we present the values of σ
x
z and out-of-plane displacement for various values of exponent p of the material power-law (p
=
0,0.5,1
,
4,10) and various thickness to side ratios (a
/
h
=
4,10,100) according to the present higher-order theory using the 5-order plate theory. Results are tabulated and compared with available references. It demonstrates the results are in excellent agreement with the meshfree results [44] and other reference values under various thickness to span ratio (a
/
h) and four and five different kinds of parameters (p).Table 6
Square B-type plate in bending with consideringσ
z
x and w use the present higher-order plate theory.
p
a
/
h
σ
z
x
(
h
/
6
)
w
(
0
)
4
10
100
4
10
100
0
Meshfree [37]
0.2208
0.2227
0.2228
0.4447
0.3711
0.3568
Present
0.2392
0.2411
0.2415
0.4655
0.3748
0.3579
0.5
Meshfree [37]
0.2546
0.2581
0.2585
0.6168
0.5238
0.5058
Present
0.2406
0.2417
0.2418
0.6409
0.5240
0.5023
1
Ref. [41]
0.2613
0.2605
0.2603
0.7628
0.6324
0.6072
CLPT [37]
0.0000
0.0000
0.0000
0.6070
0.6070
0.6070
FSDT (k
=
5
/
6) [37]
0.2458
0.2458
0.2458
0.7738
0.6337
0.6073
Ref. [39] N
=
4
0.2604
0.2594
0.2593
0.7628
0.6324
0.6072
Ref. [40]
0.2742
0.2788
0.2793
0.7416
0.6305
0.6092
Meshfree [37]
0.2745
0.2789
0.2795
0.7417
0.6305
0.6092
Present
0.2287
0.2290
0.2290
0.7766
0.6339
0.6073
4
Ref. [41]
0.2429
0.2431
0.2432
1.0934
0.8321
0.7797
CLPT [37]
0.0000
0.0000
0.0000
0.7792
0.7792
0.7792
FSDT (k
=
5
/
6) [37]
0.1877
0.1877
0.1877
1.0285
0.8191
0.7796
Ref. [39] N
=
4
0.2400
0.2398
0.2398
1.0930
0.8307
0.7797
Ref. [40]
0.2723
0.2778
0.2785
1.0391
0.8202
0.7784
Meshfree [37]
0.2696
0.2747
0.2753
1.0371
0.8199
0.7784
Present
0.1850
0.1861
0.1862
1.1034
0.8311
0.7797
10
Ref. [41]
0.2150
0.2174
0.2179
1.2232
0.8753
0.8077
CLPT [37]
0.0000
0.0000
0.0000
0.8070
0.8070
0.8070
FSDT (k
=
5
/
6) [37]
0.1234
0.1234
0.1234
1.1109
0.8556
0.8075
Ref. [39] N
=
4
0.1932
0.1944
0.1946
1.2172
0.8740
0.8077
Ref. [40]
0.2016
0.2059
0.2064
1.1780
0.8650
0.8050
Meshfree [37]
0.1995
0.2034
0.2039
1.1752
0.8645
0.8050
Present
0.1835
0.1862
0.1866
1.2372
0.8755
0.8077In Figure5, we also present the evolution of the displacement and stresses across the thickness direction according to present shear deformation theory for various values of the exponent p of a plate with side-to-thickness ratio a
/
h
=
100, using the 5-order plate theory. Compared with the meshfree method, all the values are consisted with each other except the late one. The major difference is located in internal boundary (h
=
2
h
/
5) on out-of-plane stress σ
z
x. In addition, the present solution of the stresses σ
z
x is very close to zero at upper skin and bottom skin (p
=
1,4
,
10), and the meshfree solution is away from the origin of axis. In conclusion, these results indicate that the OHSNDPT is a good choice to deal with the FGM plate problem.Figure 5
Square B-type plate subjected to sinusoidal load at the top, witha
/
h
=
100. Dimensionless displacement (w
¯) through the thickness direction according to present higher-order theory for different values of p.


### 5.1.3. Isotropic FGM Plate Compared with Laminated Plate

To demonstrate the availability of the present orthogonal higher-order shear and deformation plate theory, more cases of the modulus form are considered, linear model, exponential model, and reciprocal model (refer to [47]).The attention is focused on the influence of the three distributions of the elastic modulus on the displacement and stress fields in the plate. Young’s moduli at the upper and lower surfaces of the plate are given asE
(
0
)
=
1
G
P
a and E
(
h
)
=
10
G
P
a. Accordingly, material parameters for FGM models used in the present study are determined as follows:(a)
Liner modelE
(
z
)
=
E
0
(
p
(
z
/
h
+
1
/
2
)
+
1
), E
0
=
E
(
0
)
=
1
G
P
a, p
=
9
(b)
Exponential model,E
(
z
)
=
E
0
e
p
∗
(
z
/
h
+
1
/
2
), E
0
=
E
(
0
)
=
1
G
P
a, p
=
ln
⁡
E
h
/
E
0
=
ln
⁡
10
(c)
Reciprocal model,E
(
z
)
=
E
0
/
p
∗
(
z
/
h
+
1
/
2
)
+
1, E
0
=
E
(
0
)
=
1
G
P
a, p
=
E
(
0
)
/
E
(
h
)
-
1
=
-
0.9Based on the above three models of Young’s modulus, the variation of displacements u and w, stresses σ
x
,
τ
x
y
,
σ
z, and τ
z
x at a chosen location (
x
/
a
=
1
/
4
,
y
/
b
=
1
/
4
), as a function of the z-coordinate, are shown in Figure 6 for the plate (
h
/
a
=
0.1
).Figure 6
D-type square plate for several exponents for the powerp compared with the type E schemes.
(a)
(b)
(c)
(d)From Figure6, the following observations can be made:(I)
The magnitude of the vertical displacement (deflection) for the reciprocal model is biggest while that of the linear model is the least (Figure6(a)). This reveals that the bending rigidities of the plate are different for these three graded models. The linear model is most rigid and the reciprocal model is softest.
(II)
The in-plane stresses,σ
x
x (Figure 6(c)) and τ
x
y (Figure 6(d)), are essentially nonlinear over the thickness of the plate for all these three models of functionally graded materials. The in-plane stress concentrations in the plate are quite different for different grade models of Young’s modulus. For example, the maximum compressive stress for the reciprocal model is about 1.8 times that of the linear model (see Figure 6(c), especially 63.0921 and 34.2402 at the center point). This enables an optimal design of the plate by selecting appropriately graded models of materials.
(III)
Compared with the 3D closed-form solutions [47], the present solutions of the OHOSNDPT are well coincident with the reference data from Figure 6. Not only the in-plane stresses, but also σ
x and τ
x
y reach an identical with the 3D closed-form solutions, but also the out-of-plane stress, τ
z
x, is in very good agreement with the similar values in the literature.In short, the present approach can actually reflect the distribution situation of displacements and stresses. In addition, it is also in very good agreement with similar theories in the literature. Thus, the present method is a good idea worth trying.


### 5.1.4. Power-Law Fitted with the Straight Line

In this section, we state another scheme to illustrate the superiority of the OHOSNDPT.Because any function can be fitted with straight lines, the stiffness functionψ
(
z
) of the FGM plate with arbitrary form equivalents to a piecewise function ψ
¯
(
z
) with the function is a linear function ψ
¯
k
(
z
) in each interval. Moreover, the present orthogonal higher-order shear and nonamial plate theory also be suitable for the piecewise function ψ
¯
(
z
) in laminated FGM plates.As shown in Figure7, the variations of Young’s modulus are along the thickness of the power-law FGM plates although their values at the upper and lower surfaces of the plate are fixed and the corresponding piecewise linear functions ψ
¯
(
z
) are chosen.Figure 7
D-type square plate for several exponents for the powerp compared with the type E schemes.Constant coefficients are usually used in the each layer when the FGM plate can be seen as the laminated plates [47]. This scheme facilitates settling the FGM models, but some shortcomings are also exposed in Figure 8. Figure 8 shows the variation of the longitudinal displacements w, in-plane stresses σ
x
x and σ
x
y and out-plane stresses σ
y
z of the Tape A simple supported FGM square plate (k
=
3, p
=
10 and a
/
h
=
10).(i)
The constant coefficients of the FGM material properties in each layer are generally taken to be the value of the midpoint of the layer interval. For the class of material which the volume fraction varies steeply (p
=
10), the errors between the laminated plate models and the FGM properties obviously exist. Thus, the errors between the scheme of straight line and the analytical solution are smaller than the constant coefficients ones and it is verified from Figure 8(a).
(ii)
If the constant coefficients of the laminated plates model are used, the paramount material properties, the expression of Young’s modulus, and the density function are not continuous functions on vertical directionz. The displacement model is usually supposed as the continuous function on z, the discontinuous of stresses are obtained according to geometric equations and physics equations. Obviously, this is contrary to reality. From Figures 8(b) and 8(c), it can be observed that the scheme of straight line can match the analytical solution. The curves of the constant coefficients appear as interrupted phenomena on h
/
6, and the error enlarges with the increasing of z. On transverse stress σ
y
z (Figure 8(d)), the analytical solution and the schemes of the straight lines and the constant coefficients all occur the discontinuous situation on h
/
6, but the performance of the straight lines is better than the constant coefficients and close to the analytical solutions.Figure 8
D-type square plate for several exponents for the powerp compared with the type E schemes.
(a)
(b)
(c)
(d)From Table7, some phenomena can be obtained:(i)
Whenk is small (k
=
2,3
,
5), the scheme of the straight lines is better than the approach of the constant coefficients.
(ii)
Whenk is large (k
=
200,400), all two schemes achieve the convergence solutions.
(iii)
The two scheme coupled with OHOSNDPT are superior to the common scheme of the constant coefficients, while the familiar approach takesk
=
1000 to get the stable values.Table 7
The convergence study for the laywise number of the B-type plate using higher-order plate theory,p
=
10, and a
/
h
=
4.
Number
2
3
5
10
50
100
200
400
Analytical
Laywise (Error)
1.8480
1.6586
1.5094
1.4311
1.4038
1.4029
1.4027
1.4026
1.4026
31.76%
18.25%
7.61%
2.03%
0.09%
0.02%
0.01%
0.00%
Linear (Error)
1.1328
1.2422
1.3206
1.3770
1.4015
1.4023
1.4025
1.4026
19.24%
11.44%
5.85%
1.83%
0.08%
0.02%
0.01%
0.00%We draw a conclusion from the above example, the OHOSNDPT is a very suitable FGM plate theory.


## 5.1.1. Isotropic FGM Square Plate

In this example, an isotropic FGM square plate of type-A is considered. The plate is graded from aluminumE
m
=
70
G
P
a at the bottom to alumina E
c
=
380
G
P
a at the top. The transverse displacement, the normal stresses, and the in-plane and transverse shear stresses are presented in normalized form as(35)
w
¯
=
10
h
3
E
c
q
a
4
w
,
σ
¯
x
x
=
h
2
q
a
2
σ
x
x
.An initial convergence study was performed forσ
x
x
h
/
3 and transverse displacement w
(
0
) at the center of the plate, considering p
=
1, a
/
h
=
10, and the order K of the plate theory of 1, 3, 5, and 7. Compared with the references [37] and [39], results are presented in Tables 1 and 2. As seen in these tables, it is sufficient to use K
=
5 (we consider that K
=
5 should be used in the following computation in this section).Table 1
w
- convergence study for the bending analysis of A-type plate using higher-order plate theory, p
=
1, and a
/
h
=
10.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT(
k
=
5
/
6
) [37]
[39]
Meshfree [37]
w
-
0.5846
0.5891
0.5891
0.5891
0.5889
0.5890
0.5868Table 2
σ
-
x
x convergence study for the bending analysis of A-type plate using higher-order plate theory, p
=
1, and a
/
h
=
10.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (k
=
5
/
6) [37]
[39]
Meshfree [37]
σ
-
x
x
1.4968
1.4904
1.4908
1.4908
2.0150
1.5064
1.4917In Table3, we present results for σ
¯
x
x and transverse displacement for various exponents p of the power-law (18) considering the order of the plate theory equal 7. The considered side-to-thickness ratios (
a
/
h
) are 4, 10, and 100, which means thickness h equals 0.25, 0.1, and 0.01, respectively. Results are compared with the classical plate theory (CLPT), the first-order shear deformation theory (FSDT) with a correction factor k
=
5
/
6, and those from Zenkour’s generalized shear deformation theory [38], and those from Carrera et al. [25, 39], and Neves et al. [40].Table 3
A-type plate in bending. Theσ
z
x and w according to present higher-order plate theory.
p
a
/
h
σ
x
x
(
h
/
3
)
w
(
0
)
4
10
100
4
10
100
0
Meshfree [37]
≠0
0.5278
1.3176
13.161
0.3665
0.2942
0.2803
Present
0.5144
1.3123
13.1713
0.3800
0.2961
0.2804
0.5
Meshfree [37]
≠0
0.5860
1.4680
14.673
0.5493
0.4548
0.4365
Present
≠0
0.5724
1.4591
14.6446
0.5678
0.4538
0.4326
1
Ref. [41]
≠0
0.6221
14.969
14.969
0.7171
0.5875
0.5625
CLPT [37]
0
0.8060
2.0150
20.150
0.5623
0.5623
0.5623
FSDT (k
=
5
/
6) [37]
0
0.8060
2.0150
20.150
0.7291
0.5889
0.5625
Ref. [39] N
=
4
≠0
0.6221
1.5064
14.969
0.7171
0.5875
0.5625
Ref. [40]
≠0
0.5925
1.4945
14.969
0.6997
0.5845
0.5624
Meshfree [37]
≠0
0.5911
1.4917
14.945
0.7020
0.5868
0.5647
Present
0.5841
1.4908
14.9677
0.7318
0.5891
0.5625
4
Ref. [41]
≠0
0.4877
1.1971
11.923
1.1585
0.8821
0.8286
CLPT [37]
0
0.6420
1.6049
16.049
0.8281
0.8281
0.8281
FSDT (k
=
5
/
6) [37]
0
0.6420
1.6049
16.049
1.1125
0.8736
0.828
Ref. [39] N
=
4
≠0
0.4877
1.1971
11.923
1.1585
0.8821
0.8286
Ref. [40]
≠0
0.2723
0.2778
0.2785
1.0391
0.8202
0.7784
Meshfree [37]
≠0
0.4330
1.1588
11.737
1.1108
0.8700
0.8240
Present
0.4451
1.1787
11.9209
1.1732
0.8831
0.8287
10
Ref. [41]
≠0
0.3695
0.8965
8.9077
1.3745
1.0072
0.9361
CLPT [37]
0
0.4796
1.1990
11.990
0.9354
0.9354
0.9354
FSDT (k
=
5
/
6) [37]
0
0.4796
1.1990
11.990
1.3178
0.9966
0.9360
Ref. [39] N
=
4
≠0
0.1478
0.8965
8.9077
1.3745
1.0072
0.9361
Ref. [40]
≠0
0.3227
1.1783
11.932
1.3490
0.8750
0.8286
Meshfree [37]
≠0
0.3097
0.8462
8.6010
1.3334
0.9888
0.9227
Present
≠0
0.3238
0.8766
8.9058
1.4026
1.0094
0.9362Date shown in Table3 clearly illustrated, a subtle difference appears in the moderately thick plate between them and the present higher-order plate theory are in good agreement with those from the Refs. When the side-to-thickness ratio is 100, the values, whether it is on transverse displacements or stresses, are almost the same as the values of [50].In Figure4, we present the evolution of the displacement and stresses across the thickness direction according to present shear deformation theory for various values of the exponent p and side-to-thickness ratio a
/
h
=
4, using the 5-order plate theory. Compared with the results of the meshfree method in [44], it shows that the values of the in-plane stress σ
x
x, σ
x
y and out-plane stress σ
z
x are consistent with the reference value. Strangeness is the presence of a significant difference between the present theory (K
=
5) and the meshless solutions on the longitudinal displacement w on p
=
10. From Figure 4, we find out the relative error of displacement w can up to 9% and the maximum relative error of stresses between them is less than 3%. This discord phenomenon is due to the computational error and the instability of the numerical method. On the other hand, it demonstrates the analytical solution of the orthogonal higher-order shear deformable plate theory is necessary and feasible.Figure 4
A-type square plate subjected to sinusoidal load at the top, witha
/
h
=
4, nondimensional displacement (w
¯) and stresses through the thickness to present OHOSNDPT for different values of p.In summary, it can be concluded that the present higher-order theory provides excellent solution for the isotropic FGM plates.


## 5.1.2. Sandwich with FGM Core

In this example, we analyze the bending of a square sandwich B-type plate with thicknessh. The bottom skin is aluminum (
E
m
=
70
G
P
a
) with thickness h
b
=
0.1
h and the top skin is alumina (
E
c
=
380
G
P
a
) with thickness h
t
=
0.1
h. The core is in FGM with volume fraction of the ceramic according to (2).The transverse displacement and the normal stresses are presented in normalized form as(36)
τ
¯
x
y
=
-
h
2
q
a
2
τ
x
y
x
=
a
,
y
=
a
,
z
=
0
,
z
=
h
/
2
,
τ
¯
z
x
=
h
q
a
τ
z
x
x
=
a
,
y
=
a
,
z
=
-
h
/
N
e
/
2
.The convergence study was performed forσ
x
z and transverse displacement w
(
0
) considering p
=
4, a
/
h
=
100, and the order K of the plate theory of 1, 3, 5, and 7. Results are presented in Tables 4 and 5. We consider that the 5-order plate theory should be used in the following computation.Table 4
w convergence study for the bending analysis of B-type plate using higher-order plate theory, p
=
4, and a
/
h
=
100.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (K
=
5
/
6) [37]
Ref. [39]
Ref. [41]
w
0.7796
0.7797
0.7797
0.7797
0.7796
0.7797
0.7797Table 5
σ
x
z convergence study for the bending analysis of B-type plate using higher-order plate theory, p
=
4, and a
/
h
=
100.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT (K
=
5
/
6) [37]
Ref. [39]
Ref. [41]
σ
x
z
0.1176
0.1949
0.1862
0.1819
0.1877
0.2398
0.2432In Table6 we present the values of σ
x
z and out-of-plane displacement for various values of exponent p of the material power-law (p
=
0,0.5,1
,
4,10) and various thickness to side ratios (a
/
h
=
4,10,100) according to the present higher-order theory using the 5-order plate theory. Results are tabulated and compared with available references. It demonstrates the results are in excellent agreement with the meshfree results [44] and other reference values under various thickness to span ratio (a
/
h) and four and five different kinds of parameters (p).Table 6
Square B-type plate in bending with consideringσ
z
x and w use the present higher-order plate theory.
p
a
/
h
σ
z
x
(
h
/
6
)
w
(
0
)
4
10
100
4
10
100
0
Meshfree [37]
0.2208
0.2227
0.2228
0.4447
0.3711
0.3568
Present
0.2392
0.2411
0.2415
0.4655
0.3748
0.3579
0.5
Meshfree [37]
0.2546
0.2581
0.2585
0.6168
0.5238
0.5058
Present
0.2406
0.2417
0.2418
0.6409
0.5240
0.5023
1
Ref. [41]
0.2613
0.2605
0.2603
0.7628
0.6324
0.6072
CLPT [37]
0.0000
0.0000
0.0000
0.6070
0.6070
0.6070
FSDT (k
=
5
/
6) [37]
0.2458
0.2458
0.2458
0.7738
0.6337
0.6073
Ref. [39] N
=
4
0.2604
0.2594
0.2593
0.7628
0.6324
0.6072
Ref. [40]
0.2742
0.2788
0.2793
0.7416
0.6305
0.6092
Meshfree [37]
0.2745
0.2789
0.2795
0.7417
0.6305
0.6092
Present
0.2287
0.2290
0.2290
0.7766
0.6339
0.6073
4
Ref. [41]
0.2429
0.2431
0.2432
1.0934
0.8321
0.7797
CLPT [37]
0.0000
0.0000
0.0000
0.7792
0.7792
0.7792
FSDT (k
=
5
/
6) [37]
0.1877
0.1877
0.1877
1.0285
0.8191
0.7796
Ref. [39] N
=
4
0.2400
0.2398
0.2398
1.0930
0.8307
0.7797
Ref. [40]
0.2723
0.2778
0.2785
1.0391
0.8202
0.7784
Meshfree [37]
0.2696
0.2747
0.2753
1.0371
0.8199
0.7784
Present
0.1850
0.1861
0.1862
1.1034
0.8311
0.7797
10
Ref. [41]
0.2150
0.2174
0.2179
1.2232
0.8753
0.8077
CLPT [37]
0.0000
0.0000
0.0000
0.8070
0.8070
0.8070
FSDT (k
=
5
/
6) [37]
0.1234
0.1234
0.1234
1.1109
0.8556
0.8075
Ref. [39] N
=
4
0.1932
0.1944
0.1946
1.2172
0.8740
0.8077
Ref. [40]
0.2016
0.2059
0.2064
1.1780
0.8650
0.8050
Meshfree [37]
0.1995
0.2034
0.2039
1.1752
0.8645
0.8050
Present
0.1835
0.1862
0.1866
1.2372
0.8755
0.8077In Figure5, we also present the evolution of the displacement and stresses across the thickness direction according to present shear deformation theory for various values of the exponent p of a plate with side-to-thickness ratio a
/
h
=
100, using the 5-order plate theory. Compared with the meshfree method, all the values are consisted with each other except the late one. The major difference is located in internal boundary (h
=
2
h
/
5) on out-of-plane stress σ
z
x. In addition, the present solution of the stresses σ
z
x is very close to zero at upper skin and bottom skin (p
=
1,4
,
10), and the meshfree solution is away from the origin of axis. In conclusion, these results indicate that the OHSNDPT is a good choice to deal with the FGM plate problem.Figure 5
Square B-type plate subjected to sinusoidal load at the top, witha
/
h
=
100. Dimensionless displacement (w
¯) through the thickness direction according to present higher-order theory for different values of p.


## 5.1.3. Isotropic FGM Plate Compared with Laminated Plate

To demonstrate the availability of the present orthogonal higher-order shear and deformation plate theory, more cases of the modulus form are considered, linear model, exponential model, and reciprocal model (refer to [47]).The attention is focused on the influence of the three distributions of the elastic modulus on the displacement and stress fields in the plate. Young’s moduli at the upper and lower surfaces of the plate are given asE
(
0
)
=
1
G
P
a and E
(
h
)
=
10
G
P
a. Accordingly, material parameters for FGM models used in the present study are determined as follows:(a)
Liner modelE
(
z
)
=
E
0
(
p
(
z
/
h
+
1
/
2
)
+
1
), E
0
=
E
(
0
)
=
1
G
P
a, p
=
9
(b)
Exponential model,E
(
z
)
=
E
0
e
p
∗
(
z
/
h
+
1
/
2
), E
0
=
E
(
0
)
=
1
G
P
a, p
=
ln
⁡
E
h
/
E
0
=
ln
⁡
10
(c)
Reciprocal model,E
(
z
)
=
E
0
/
p
∗
(
z
/
h
+
1
/
2
)
+
1, E
0
=
E
(
0
)
=
1
G
P
a, p
=
E
(
0
)
/
E
(
h
)
-
1
=
-
0.9Based on the above three models of Young’s modulus, the variation of displacements u and w, stresses σ
x
,
τ
x
y
,
σ
z, and τ
z
x at a chosen location (
x
/
a
=
1
/
4
,
y
/
b
=
1
/
4
), as a function of the z-coordinate, are shown in Figure 6 for the plate (
h
/
a
=
0.1
).Figure 6
D-type square plate for several exponents for the powerp compared with the type E schemes.
(a)
(b)
(c)
(d)From Figure6, the following observations can be made:(I)
The magnitude of the vertical displacement (deflection) for the reciprocal model is biggest while that of the linear model is the least (Figure6(a)). This reveals that the bending rigidities of the plate are different for these three graded models. The linear model is most rigid and the reciprocal model is softest.
(II)
The in-plane stresses,σ
x
x (Figure 6(c)) and τ
x
y (Figure 6(d)), are essentially nonlinear over the thickness of the plate for all these three models of functionally graded materials. The in-plane stress concentrations in the plate are quite different for different grade models of Young’s modulus. For example, the maximum compressive stress for the reciprocal model is about 1.8 times that of the linear model (see Figure 6(c), especially 63.0921 and 34.2402 at the center point). This enables an optimal design of the plate by selecting appropriately graded models of materials.
(III)
Compared with the 3D closed-form solutions [47], the present solutions of the OHOSNDPT are well coincident with the reference data from Figure 6. Not only the in-plane stresses, but also σ
x and τ
x
y reach an identical with the 3D closed-form solutions, but also the out-of-plane stress, τ
z
x, is in very good agreement with the similar values in the literature.In short, the present approach can actually reflect the distribution situation of displacements and stresses. In addition, it is also in very good agreement with similar theories in the literature. Thus, the present method is a good idea worth trying.


## 5.1.4. Power-Law Fitted with the Straight Line

In this section, we state another scheme to illustrate the superiority of the OHOSNDPT.Because any function can be fitted with straight lines, the stiffness functionψ
(
z
) of the FGM plate with arbitrary form equivalents to a piecewise function ψ
¯
(
z
) with the function is a linear function ψ
¯
k
(
z
) in each interval. Moreover, the present orthogonal higher-order shear and nonamial plate theory also be suitable for the piecewise function ψ
¯
(
z
) in laminated FGM plates.As shown in Figure7, the variations of Young’s modulus are along the thickness of the power-law FGM plates although their values at the upper and lower surfaces of the plate are fixed and the corresponding piecewise linear functions ψ
¯
(
z
) are chosen.Figure 7
D-type square plate for several exponents for the powerp compared with the type E schemes.Constant coefficients are usually used in the each layer when the FGM plate can be seen as the laminated plates [47]. This scheme facilitates settling the FGM models, but some shortcomings are also exposed in Figure 8. Figure 8 shows the variation of the longitudinal displacements w, in-plane stresses σ
x
x and σ
x
y and out-plane stresses σ
y
z of the Tape A simple supported FGM square plate (k
=
3, p
=
10 and a
/
h
=
10).(i)
The constant coefficients of the FGM material properties in each layer are generally taken to be the value of the midpoint of the layer interval. For the class of material which the volume fraction varies steeply (p
=
10), the errors between the laminated plate models and the FGM properties obviously exist. Thus, the errors between the scheme of straight line and the analytical solution are smaller than the constant coefficients ones and it is verified from Figure 8(a).
(ii)
If the constant coefficients of the laminated plates model are used, the paramount material properties, the expression of Young’s modulus, and the density function are not continuous functions on vertical directionz. The displacement model is usually supposed as the continuous function on z, the discontinuous of stresses are obtained according to geometric equations and physics equations. Obviously, this is contrary to reality. From Figures 8(b) and 8(c), it can be observed that the scheme of straight line can match the analytical solution. The curves of the constant coefficients appear as interrupted phenomena on h
/
6, and the error enlarges with the increasing of z. On transverse stress σ
y
z (Figure 8(d)), the analytical solution and the schemes of the straight lines and the constant coefficients all occur the discontinuous situation on h
/
6, but the performance of the straight lines is better than the constant coefficients and close to the analytical solutions.Figure 8
D-type square plate for several exponents for the powerp compared with the type E schemes.
(a)
(b)
(c)
(d)From Table7, some phenomena can be obtained:(i)
Whenk is small (k
=
2,3
,
5), the scheme of the straight lines is better than the approach of the constant coefficients.
(ii)
Whenk is large (k
=
200,400), all two schemes achieve the convergence solutions.
(iii)
The two scheme coupled with OHOSNDPT are superior to the common scheme of the constant coefficients, while the familiar approach takesk
=
1000 to get the stable values.Table 7
The convergence study for the laywise number of the B-type plate using higher-order plate theory,p
=
10, and a
/
h
=
4.
Number
2
3
5
10
50
100
200
400
Analytical
Laywise (Error)
1.8480
1.6586
1.5094
1.4311
1.4038
1.4029
1.4027
1.4026
1.4026
31.76%
18.25%
7.61%
2.03%
0.09%
0.02%
0.01%
0.00%
Linear (Error)
1.1328
1.2422
1.3206
1.3770
1.4015
1.4023
1.4025
1.4026
19.24%
11.44%
5.85%
1.83%
0.08%
0.02%
0.01%
0.00%We draw a conclusion from the above example, the OHOSNDPT is a very suitable FGM plate theory.


## 5.2. Free Vibrations

In this example, we study the free vibration of a simply supported isotropic FGM square plate (a
=
b
=
1) of type A. The plate is graded from aluminum (bottom) to zirconia/alumina (top). The values for metal and ceramics in the FGM plate are listed in Table 8.Table 8
Material properties of the used FGM plate.
Material
Properties
E
(
G
P
a
)
ν
ρ
(
k
g
/
m
3
)
Aluminum (A
l)
70
0.3
2702
Alumina (A
l
2
O
3)
380
0.3
3800
Zirconia (Z
r
O
2)
200
0.3
5700We consider the Mori–Tanaka homogenization scheme (36), as in Vel and Batra [51] (here considered to be the exact solution), and as in Neves et al. [40] and Qian et al. [44].The frequencyw has been nondimensionalized as follows:(37)
w
¯
=
w
h
ρ
m
E
m
.


### 5.2.1. Isotropic FGM Plate

In Table9 we present the results obtained with the theories considered and different values of p for a side-to-thickness ratio a
/
h
=
5. Results show that the convergence values are obtained with K
=
3. It is noteworthy that the present solutions are almost the same as values of [48] for p
=
1 regardless of the ratio length-thickness.Table 9
Fundamental frequency of a square simply supported A-type square plate (Al/ZrO2) using the present higher-order theory.
Mode
p
=
1
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
p
=
2
p
=
3
p
=
5
Exact 3D [48]
0.0153
0.0596
0.2192
0.2197
0.2211
0.2225
HOSNT-15 [49]
0.0154
0.0596
0.2191
0.2196
0.2211
0.2225
Ref. [37]
0.0153
0.0596
0.2193
0.2200
0.2215
0.2230
Ref. [44]
0.0149
0.0584
0.2152
0.2153
0.2172
0.2194
Ref. [42]
0.0158
0.0619
0.2276
0.2264
0.2276
0.2291
Present (K
=
1)
0.0159
0.0621
0.2307
0.2295
0.2309
0.2326
Present (K
=
3)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2272
Present (K
=
5)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2271
Present (K
=
7)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2271To further compare the exact HSDT solution and the present OHSNDPT, we introduce more reliable data. In Table10, based on the present closed-form solutions, exact HSDT solution with refined plate theory and finite element method, numerical results have been performed for simple supported Al/Al2O3 square plates when p
=
1. A well-known commercially available FEM package was used to investigate 3D free vibration of FGM square plates in [48].Table 10
The natural frequency of A-type square simply supported plate (A
l
/
A
l
2
O
3) using the higher-order theory.
a
/
h
Model
1
(
1,1
)
2
(
1,2
)
3
(
2,1
)
4
(
2,2
)
5
(
1,3
)
6
(
3,1
)
7
(
2,3
)
8
(
3,2
)
9
(
1,4
)
1
0
(
4,1
)
20
Present
359.90
889.22
889.22
1406.6
1745.3
1745.3
2244.4
2244.4
2894.3
2894.3
Ref. [43]
359.97
889.39
889.39
1406.9
1745.7
1745.7
2244.9
2244.9
—
—
FEM [43]
357.37
883.58
883.58
1398.6
1736.2
1736.2
2234.1
2234.1
—
—
Error
-
0.71
%
-
0.64
%
-
0.64
%
-
0.57
%
-
0.52
%
-
0.52
%
-
0.46
%
-
0.46
%
10
Present
703.30
1685.4
1685.4
2596.2
3170.5
3170.5
3990.0
3990.0
5016.1
5016.1
Ref. [43]
703.45
1685.7
1685.7
2596.6
3170.7
3170.7
3989.9
—
—
—
FEM [43]
699.30
1679.7
1679.7
2592.5
3169.6
3169.6
3994.8
—
—
—
Error
-
0.57
%
-
0.34
%
-
0.34
%
-
0.14
%
-
0
.
0
3
%
-
0.03
%
0.12%
5
Present
1298.1
2871.7
2871.7
4186.8
4969.4
4969.4
6041.0
6041.0
7327.4
7327.4
Ref. [43]
1298.3
2870.5
2870.5
4181.8
4960.5
4960.5
—
—
—
—
FEM [43]
1296.3
2883.0
2883.0
4216.1
5010.6
5010.6
—
—
—
—
Error
-
0.14
%
0.39%
0.39%
0.69%
0.82%
0.82%The length of square plates is 1 m. Three different thicknesses 0.05 m (corresponding to thin plates), 0.1 m and 0.2 m (corresponding to moderately thick plates) have been used. The calculations of the present scheme are obtained for first ten natural frequencies. “—” values in Table10 indicate no calculation in corresponding reference. The percentage difference given in Table 10 is defined as follows: E
r
r
o
r
(
%
)
=
[
F
E
M
-
(
s
o
l
u
t
i
o
n
)
]
/
F
E
M.An excellent agreement is observed between the present plate theory, exact HSDT solution and the FEM. It is seen that the exact results of proposed method are close to the FEM analysis results at lower and higher frequencies. Frequencies rise with an increase in the thickness of plate. This phenomena originates from the increasing the rigidity of plate. Frequencies decrease when less restraining boundary is used at the edges of square plates. This is due to the fact that higher constraints at the edges increase the flexural rigidity of the plate, leading to a higher frequency response. Compared with the exact HSDT solution, the present closed-form solution is closer to the FEM solution.


### 5.2.2. Isotropic FGM Plate Compared with Laminated Plate

In this subsection the higher-order plate theories are used for the free vibration analysis of simply supported functionally of three models of type D with side-to-thickness ratioa
/
h
=
10. All the material properties is the same as that in Section 5.1.3 and the density function is supposed with ρ
(
z
)
=
ρ
0
ψ
(
z
), that is, ρ
(
h
/
2
)
=
10
∗
ρ
(
-
h
/
2
).In Table11, nondimensional first 10 frequencies of plate are computed for three different parameter a
/
h
=
20,10,5 in the Tape D. Some views can be noticeable:(i)
The order of the first ten frequency of the three model all comply with the order of isotropic square plate (m
2
+
n
2). It shows the variety of the properties of material on z has nothing to do with the order of the frequencies.
(ii)
Nondimensional results are hardly affected bya
/
h.
(iii)
Unlike static problem, the longitudinal displacementw for the reciprocal model is about 1.8 times that of the linear model. The frequencies of them are very similar to each other for fixed (
m
,
n
). It is because the ratio of the upper material and the bottom material of the density function is set as 10, which is consistent with the ratio of elastic modulus.Table 11
First ten frequencies of type D square plate with variousa
/
h using the higher-order theory.
Source
1
(
1,1
)
2
(
1,2
)
3
(
2,1
)
4
(
2,2
)
5
(
1,3
)
6
(
3,1
)
7
(
2,3
)
8
(
3,2
)
9
(
1,4
)
1
0
(
4,1
)
Linear model
a
/
h
=
20
0.8481
2.0986
2.0986
3.3243
4.1284
4.1284
5.3158
5.3158
6.8659
6.8659
a
/
h
=
10
0.8311
2.0012
2.0012
3.0947
3.7874
3.7874
4.7796
4.7796
6.0269
6.0269
a
/
h
=
5
0.7737
1.7284
1.7284
2.5333
3.0136
3.0136
3.6719
3.6719
4.4619
4.4619
Exponent model
a
/
h
=
20
0.8504
2.1025
2.1025
3.3280
4.1310
4.1310
5.3155
5.3155
6.8596
6.8596
a
/
h
=
10
0.8320
1.9982
1.9982
3.0833
3.7688
3.7688
4.7485
4.7485
5.9769
5.9769
a
/
h
=
5
0.7708
1.7121
1.7121
2.5011
2.9708
2.9708
3.6139
3.6139
4.3853
4.3853
Reciprocal model
a
/
h
=
20
0.9214
2.2722
2.2722
3.5877
4.4464
4.4464
5.7559
5.7559
7.3460
7.3460
a
/
h
=
10
0.8969
2.1357
2.1357
3.2725
3.9843
3.9843
4.9939
4.9939
6.2496
6.2496
a
/
h
=
5
0.8181
1.7836
1.7836
2.5767
3.0444
3.0444
3.6811
3.6811
4.4410
4.4410


## 5.2.1. Isotropic FGM Plate

In Table9 we present the results obtained with the theories considered and different values of p for a side-to-thickness ratio a
/
h
=
5. Results show that the convergence values are obtained with K
=
3. It is noteworthy that the present solutions are almost the same as values of [48] for p
=
1 regardless of the ratio length-thickness.Table 9
Fundamental frequency of a square simply supported A-type square plate (Al/ZrO2) using the present higher-order theory.
Mode
p
=
1
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
p
=
2
p
=
3
p
=
5
Exact 3D [48]
0.0153
0.0596
0.2192
0.2197
0.2211
0.2225
HOSNT-15 [49]
0.0154
0.0596
0.2191
0.2196
0.2211
0.2225
Ref. [37]
0.0153
0.0596
0.2193
0.2200
0.2215
0.2230
Ref. [44]
0.0149
0.0584
0.2152
0.2153
0.2172
0.2194
Ref. [42]
0.0158
0.0619
0.2276
0.2264
0.2276
0.2291
Present (K
=
1)
0.0159
0.0621
0.2307
0.2295
0.2309
0.2326
Present (K
=
3)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2272
Present (K
=
5)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2271
Present (K
=
7)
0.0158
0.0619
0.2277
0.2256
0.2262
0.2271To further compare the exact HSDT solution and the present OHSNDPT, we introduce more reliable data. In Table10, based on the present closed-form solutions, exact HSDT solution with refined plate theory and finite element method, numerical results have been performed for simple supported Al/Al2O3 square plates when p
=
1. A well-known commercially available FEM package was used to investigate 3D free vibration of FGM square plates in [48].Table 10
The natural frequency of A-type square simply supported plate (A
l
/
A
l
2
O
3) using the higher-order theory.
a
/
h
Model
1
(
1,1
)
2
(
1,2
)
3
(
2,1
)
4
(
2,2
)
5
(
1,3
)
6
(
3,1
)
7
(
2,3
)
8
(
3,2
)
9
(
1,4
)
1
0
(
4,1
)
20
Present
359.90
889.22
889.22
1406.6
1745.3
1745.3
2244.4
2244.4
2894.3
2894.3
Ref. [43]
359.97
889.39
889.39
1406.9
1745.7
1745.7
2244.9
2244.9
—
—
FEM [43]
357.37
883.58
883.58
1398.6
1736.2
1736.2
2234.1
2234.1
—
—
Error
-
0.71
%
-
0.64
%
-
0.64
%
-
0.57
%
-
0.52
%
-
0.52
%
-
0.46
%
-
0.46
%
10
Present
703.30
1685.4
1685.4
2596.2
3170.5
3170.5
3990.0
3990.0
5016.1
5016.1
Ref. [43]
703.45
1685.7
1685.7
2596.6
3170.7
3170.7
3989.9
—
—
—
FEM [43]
699.30
1679.7
1679.7
2592.5
3169.6
3169.6
3994.8
—
—
—
Error
-
0.57
%
-
0.34
%
-
0.34
%
-
0.14
%
-
0
.
0
3
%
-
0.03
%
0.12%
5
Present
1298.1
2871.7
2871.7
4186.8
4969.4
4969.4
6041.0
6041.0
7327.4
7327.4
Ref. [43]
1298.3
2870.5
2870.5
4181.8
4960.5
4960.5
—
—
—
—
FEM [43]
1296.3
2883.0
2883.0
4216.1
5010.6
5010.6
—
—
—
—
Error
-
0.14
%
0.39%
0.39%
0.69%
0.82%
0.82%The length of square plates is 1 m. Three different thicknesses 0.05 m (corresponding to thin plates), 0.1 m and 0.2 m (corresponding to moderately thick plates) have been used. The calculations of the present scheme are obtained for first ten natural frequencies. “—” values in Table10 indicate no calculation in corresponding reference. The percentage difference given in Table 10 is defined as follows: E
r
r
o
r
(
%
)
=
[
F
E
M
-
(
s
o
l
u
t
i
o
n
)
]
/
F
E
M.An excellent agreement is observed between the present plate theory, exact HSDT solution and the FEM. It is seen that the exact results of proposed method are close to the FEM analysis results at lower and higher frequencies. Frequencies rise with an increase in the thickness of plate. This phenomena originates from the increasing the rigidity of plate. Frequencies decrease when less restraining boundary is used at the edges of square plates. This is due to the fact that higher constraints at the edges increase the flexural rigidity of the plate, leading to a higher frequency response. Compared with the exact HSDT solution, the present closed-form solution is closer to the FEM solution.


## 5.2.2. Isotropic FGM Plate Compared with Laminated Plate

In this subsection the higher-order plate theories are used for the free vibration analysis of simply supported functionally of three models of type D with side-to-thickness ratioa
/
h
=
10. All the material properties is the same as that in Section 5.1.3 and the density function is supposed with ρ
(
z
)
=
ρ
0
ψ
(
z
), that is, ρ
(
h
/
2
)
=
10
∗
ρ
(
-
h
/
2
).In Table11, nondimensional first 10 frequencies of plate are computed for three different parameter a
/
h
=
20,10,5 in the Tape D. Some views can be noticeable:(i)
The order of the first ten frequency of the three model all comply with the order of isotropic square plate (m
2
+
n
2). It shows the variety of the properties of material on z has nothing to do with the order of the frequencies.
(ii)
Nondimensional results are hardly affected bya
/
h.
(iii)
Unlike static problem, the longitudinal displacementw for the reciprocal model is about 1.8 times that of the linear model. The frequencies of them are very similar to each other for fixed (
m
,
n
). It is because the ratio of the upper material and the bottom material of the density function is set as 10, which is consistent with the ratio of elastic modulus.Table 11
First ten frequencies of type D square plate with variousa
/
h using the higher-order theory.
Source
1
(
1,1
)
2
(
1,2
)
3
(
2,1
)
4
(
2,2
)
5
(
1,3
)
6
(
3,1
)
7
(
2,3
)
8
(
3,2
)
9
(
1,4
)
1
0
(
4,1
)
Linear model
a
/
h
=
20
0.8481
2.0986
2.0986
3.3243
4.1284
4.1284
5.3158
5.3158
6.8659
6.8659
a
/
h
=
10
0.8311
2.0012
2.0012
3.0947
3.7874
3.7874
4.7796
4.7796
6.0269
6.0269
a
/
h
=
5
0.7737
1.7284
1.7284
2.5333
3.0136
3.0136
3.6719
3.6719
4.4619
4.4619
Exponent model
a
/
h
=
20
0.8504
2.1025
2.1025
3.3280
4.1310
4.1310
5.3155
5.3155
6.8596
6.8596
a
/
h
=
10
0.8320
1.9982
1.9982
3.0833
3.7688
3.7688
4.7485
4.7485
5.9769
5.9769
a
/
h
=
5
0.7708
1.7121
1.7121
2.5011
2.9708
2.9708
3.6139
3.6139
4.3853
4.3853
Reciprocal model
a
/
h
=
20
0.9214
2.2722
2.2722
3.5877
4.4464
4.4464
5.7559
5.7559
7.3460
7.3460
a
/
h
=
10
0.8969
2.1357
2.1357
3.2725
3.9843
3.9843
4.9939
4.9939
6.2496
6.2496
a
/
h
=
5
0.8181
1.7836
1.7836
2.5767
3.0444
3.0444
3.6811
3.6811
4.4410
4.4410


## 5.3. Buckling

### 5.3.1. Sandwich Plate with FGM Skins

In the next examples the higher-order plate theories are used for the buckling analysis of simply supported functionally graded sandwich square plates (a
=
b
=
1) of type C with side-to-thickness ratio a
/
h
=
10. The uni-and biaxial critical buckling loads are analsed.The material properties areE
m
=
70
E
0 (aluminum) for the metal and E
c
=
380
E
0 (alumina) for the ceramic being E
0
=
1
G
P
a. The law-of-mixtures ((19)–(21)) was used for the computation of Young’s modulus for each layer. The nondimensional parameter used is(38)
P
¯
=
P
a
2
100
h
3
E
0
.An initial convergence study with the higher-order theory was conducted for each buckling load type considering the orderK of the plate theory is 1, 3, 5, and 7. The uniaxial case is presented in Table 12 for the 2
-
2
-
1 sandwich with p
=
5 and the biaxial case is presented in Table 13 for the 1
-
2
-
1 sandwich with p
=
1. Further results are obtained by considering K equal to 5, which seems acceptable by the convergence study.Table 12
Convergence study for the uniaxial buckling load of a simply supported 2-2-1 sandwich square plate with FGM skins andp
=
5 case using the higher-order theory.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT [37]
TSDPT [45]
SSDPT [46]
P
-
4.11230
4.10740
4.10605
4.10449
4.09285
4.11209
4.11280Table 13
Convergence study for the uniaxial buckling load of a simply supported 1-2-1 sandwich square plate with FGM skins andp
=
1 case using the higher-order theory.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT [37]
TSDPT [45]
SSDPT [46]
P
-
3.76605
3.75329
3.75304
3.75290
3.74182
3.75328
3.75314The critical buckling loads obtained from the present approach are tabulated in Tables14 and 15 for various power-law exponents p and thickness ratios. Both tables include results obtained from classical plate theory (CLPT), first-order shear deformation plate theory (FSDPT, the shear correction factor is taken as 5
/
6), Reddy’s higher-order shear deformation plate theory (TSDPT) [10], and Zenkour’s sinusoidal shear deformation plate theory (SSDPT) [29]. Table 14 refers to the uniaxial buckling load and Table 15 refers to the biaxial buckling load.Table 14
Uniaxial buckling load of simply supported plate of C-type using the 5-order theory.
p
Theory
P
-
1-0-1
2-1-2
2-1-1
1-1-1
2-2-1
1-2-1
0
CLPT [37]
13.73791
13.73791
13.73791
13.73791
13.73791
13.73791
FSDPT [37]
13.00449
13.00449
13.00449
13.00449
13.00449
13.00449
TSDPT [45]
13.00495
13.00495
13.00495
13.00495
13.00495
13.00495
SSDPT [46]
13.00606
13.00606
13.00606
13.00606
13.00606
13.00606
Meshfree [37]
13.95287
13.95287
13.95287
13.95287
13.95287
13.95287
13.00516
13.00516
13.00516
13.00516
13.00516
13.00516
0.5
CLPT [37]
7.65398
8.25597
8.56223
8.78063
9.18254
9.61525
FSDPT [37]
7.33732
7.91320
8.20015
8.41034
8.78673
9.19517
TSDPT [45]
7.36437
7.94084
8.22470
8.43645
8.80997
9.21681
SSDPT [46]
7.36568
7.94195
8.22538
8.43712
8.81037
9.21670
Meshfree [37]
7.16207
7.71627
7.98956
8.22133
8.55172
8.94190
7.36361
7.94047
8.22048
8.43637
8.80980
9.21651
1
CLPT [37]
5.33248
6.02733
6.40391
6.68150
7.19663
7.78406
FSDPT [37]
5.14236
5.81379
6.17020
6.43892
6.92571
7.48365
TSDPT [45]
5.16713
5.84006
6.19394
6.46474
6.94944
7.50656
SSDPT [46]
5.16846
5.84119
6.19461
6.46539
6.94980
7.50629
Meshfree [37]
5.06137
5.71135
6.05467
6.31500
6.78405
7.31995
5.16489
5.83877
6.19209
6.48150
6.94876
7.50579
5
CLPT [37]
2.73080
3.10704
3.48418
3.65732
4.21238
4.85717
FSDPT [37]
2.63842
3.02252
3.38538
3.55958
4.09285
4.71475
TSDPT [45]
2.65821
3.04257
3.40351
3.57956
4.11209
4.73469
SSDPT [46]
2.66006
3.04406
3.40449
3.58063
4.11288
4.73488
Meshfree [37]
2.63652
3.00791
3.36225
3.53005
4.05070
4.64701
2.64204
3.02902
3.38599
3.57260
4.10499
4.73116
10
CLPT [37]
2.56985
2.80340
3.16427
3.25924
3.79238
4.38221
FSDPT [37]
2.46904
2.72626
3.07428
3.17521
3.68890
4.26040
TSDPT [45]
2.48727
2.74632
3.09190
3.19471
3.70752
4.27991
SSDPT [46]
2.48928
2.74844
3.13443
3.19456
3.14574
4.38175
Meshfree [37]
2.47216
2.72046
3.06067
3.15761
3.66166
4.20550
2.47035
2.72423
3.06339
3.18364
3.09531
4.27493Table 15
Biaxial buckling load of simply supported plate of C-type using the 5-order theory.
p
Theory
P
-
Source
1-0-1
2-1-2
2-1-1
1-1-1
2-2-1
1-2-1
0
CLPT [37]
6.86896
6.86896
6.86896
6.86896
6.86896
6.86896
FSDPT [37]
6.50224
6.50224
6.50224
6.50224
6.50224
6.50224
TSDPT [45]
6.50248
6.50248
6.50248
6.50248
6.50248
6.50248
SSDPT [46]
6.50303
6.50303
6.50303
6.50303
6.50303
6.50303
Meshfree [37]
6.47643
6.47643
6.47643
6.47643
6.47643
6.47643
6.50258
6.50258
6.50258
6.50258
6.50258
6.50258
0.5
CLPT [37]
3.82699
4.12798
4.28112
4.39032
4.59127
4.80762
FSDPT [37]
3.66866
3.95660
4.10007
4.20517
4.39336
4.59758
TSDPT [45]
3.68219
3.97042
4.11235
4.21823
4.40499
4.60841
SSDPT [46]
3.68284
3.97097
4.11269
4.21856
4.40519
4.60835
Meshfree [37]
3.58104
3.85813
3.99478
4.09639
4.27586
4.47095
3.68180
3.97024
4.11204
4.21819
4.40490
4.60826
1
CLPT [37]
2.66624
3.01366
3.20195
3.34075
3.59831
3.89203
FSDPT [37]
2.57118
2.90690
3.08510
3.21946
3.46286
3.74182
TSDPT [45]
2.58357
2.92003
3.09697
3.23237
3.47472
3.75328
SSDPT [46]
2.58423
2.92060
3.09731
3.23270
3.47490
3.75314
Meshfree [37]
2.53069
2.85568
3.02733
3.15750
3.39202
3.65998
2.58245
2.91939
3.09605
3.23212
3.47438
3.75290
5
CLPT [37]
1.36540
1.55352
1.74209
1.82866
2.10619
2.42859
FSDPT [37]
1.31921
1.51126
1.69269
1.77979
2.04642
2.35737
TSDPT [45]
1.32910
1.52129
1.70176
1.78978
2.05605
2.36734
SSDPT [46]
1.32910
1.52203
1.70224
1.79032
2.05644
2.36744
Meshfree [37]
1.31826
1.50395
1.68128
1.76502
2.02535
2.32351
1.32102
1.51451
1.69299
1.78630
2.05225
2.36558
10
CLPT [37]
1.28493
1.40170
1.58214
1.62962
1.89619
2.19111
FSDPT [37]
1.23452
1.36313
1.53714
1.58760
1.84445
2.13020
TSDPT [45]
1.24363
1.37316
1.54595
1.59736
1.85376
2.13995
SSDPT [46]
1.24475
1.37422
1.56721
1.59728
1.57287
2.19087
Meshfree [37]
1.23608
1.36023
1.53034
1.57880
1.83083
2.10275
1.23518
1.36211
1.53170
1.59182
1.84765
2.13746A good agreement between the present solution and references considered; especially [10, 29] is obtained. This allows us to conclude that the present higher-order plate theory is good for the modeling of simply supported sandwich FGM plates with OHOSNDPT as a good formulation. This study also lead us to conclude that the thickness stretching effect has a strong influence on the buckling analysis of sandwich FGM plates as gives higher fundamental buckling loads.The isotropic fully ceramic plate (first line on Tables14 and 15) has the higher fundamental buckling loads. As the core thickness to the total thickness of the plate ratio (
(
h
2
-
h
1
)
/
h
) increases, the buckling loads increase as well. Considering each column of both tables we may conclude that the critical buckling loads decrease as the power-law exponent p increases. By comparing Tables 14 and 15 we also conclude that the biaxial buckling load of simply supported sandwich square plate with FGM skins is half the uniaxial one for the same plate.


### 5.3.2. Isotropic FGM Plate Compared with Laminated Plate

In this subsection, the higher-order plate theories are used for the buckling analysis of simply supported functionally of two models of type D and E with side-to-thickness ratioa
/
h
=
10. All the material properties are taken as in Section 5.1.3.The convergence study with the higher-order theory was also conducted for each buckling load type considering the orderK of the plate theory from 1 to 7. The uniaxial case is presented in Table 16 for the Linear model with a
/
h
=
20,10,5. Further results are obtained by considering K equals 3, which seems acceptable by the convergence study.Table 16
Convergence study for the uniaxial buckling load of a simply supported square plate with linear model using the higher-order theory.
Scheme
K
=
1
K
=
2
K
=
3
K
=
4
K
=
5
K
=
6
K
=
7
a
/
h
=
20
30.6144
30.5958
30.5585
30.5584
30.5584
30.5584
30.5584
a
/
h
=
10
14.9028
14.8679
14.7979
14.7978
14.7978
14.7977
14.7977
a
/
h
=
5
6.7392
6.6849
6.5739
6.5739
6.5738
6.5738
6.5738In Table17, nondimensional uniaxial buckling load and biaxial buckling load of plate are computed for various a
/
h of three different models.Table 17
Bi and iniaxial buckling load of simply supported plate of type D using the higher-order theory.
Model
Reciprocal
Exponent
Linear
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
Biaxial
15.2792
7.3989
3.2869
10.9162
5.2696
2.3170
8.3944
4.0156
1.7136
Uniaxial
30.5585
14.7979
6.5739
21.8325
10.5393
4.6339
16.7888
8.0313
3.4273It is well known that as the plate thickness increases, the shear effect which leads to reduction of the critical buckling stress increases. It can be seen from Tables16 and 17 that although increasing the plate thickness decreases P
¯
x in all models, it can be concluded that for compressive biaxial and uniaxial loadings consideration of the prebuckling deformations in the stability equations increases the shear effect.Another interesting phenomenon appears between the static problem and the buckling situation. Whena
/
h
=
10, that is, the FGM is a thin plate, the rate of the buckling stresses between the reciprocal model and the linear model (the value is 1.8425) is almost the same as the proportion of the displacement at the center point (the more accurate value 1.8426 is than the number in Section 5.1.4) at the same K. When a
/
h
=
5, that is, the FGM is a moderately thick plate, the former ratio is 1.9181 and the latter one becomes 1.9198 (4.6322/2.4128). In summary, it not only illustrate that the present OHOSNDPT is high accuracy both the static problem and the buckling situation, but also reflect the stability for the thin plate and moderately thick plate.


## 5.3.1. Sandwich Plate with FGM Skins

In the next examples the higher-order plate theories are used for the buckling analysis of simply supported functionally graded sandwich square plates (a
=
b
=
1) of type C with side-to-thickness ratio a
/
h
=
10. The uni-and biaxial critical buckling loads are analsed.The material properties areE
m
=
70
E
0 (aluminum) for the metal and E
c
=
380
E
0 (alumina) for the ceramic being E
0
=
1
G
P
a. The law-of-mixtures ((19)–(21)) was used for the computation of Young’s modulus for each layer. The nondimensional parameter used is(38)
P
¯
=
P
a
2
100
h
3
E
0
.An initial convergence study with the higher-order theory was conducted for each buckling load type considering the orderK of the plate theory is 1, 3, 5, and 7. The uniaxial case is presented in Table 12 for the 2
-
2
-
1 sandwich with p
=
5 and the biaxial case is presented in Table 13 for the 1
-
2
-
1 sandwich with p
=
1. Further results are obtained by considering K equal to 5, which seems acceptable by the convergence study.Table 12
Convergence study for the uniaxial buckling load of a simply supported 2-2-1 sandwich square plate with FGM skins andp
=
5 case using the higher-order theory.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT [37]
TSDPT [45]
SSDPT [46]
P
-
4.11230
4.10740
4.10605
4.10449
4.09285
4.11209
4.11280Table 13
Convergence study for the uniaxial buckling load of a simply supported 1-2-1 sandwich square plate with FGM skins andp
=
1 case using the higher-order theory.
Scheme
K
=
1
K
=
3
K
=
5
K
=
7
FSDT [37]
TSDPT [45]
SSDPT [46]
P
-
3.76605
3.75329
3.75304
3.75290
3.74182
3.75328
3.75314The critical buckling loads obtained from the present approach are tabulated in Tables14 and 15 for various power-law exponents p and thickness ratios. Both tables include results obtained from classical plate theory (CLPT), first-order shear deformation plate theory (FSDPT, the shear correction factor is taken as 5
/
6), Reddy’s higher-order shear deformation plate theory (TSDPT) [10], and Zenkour’s sinusoidal shear deformation plate theory (SSDPT) [29]. Table 14 refers to the uniaxial buckling load and Table 15 refers to the biaxial buckling load.Table 14
Uniaxial buckling load of simply supported plate of C-type using the 5-order theory.
p
Theory
P
-
1-0-1
2-1-2
2-1-1
1-1-1
2-2-1
1-2-1
0
CLPT [37]
13.73791
13.73791
13.73791
13.73791
13.73791
13.73791
FSDPT [37]
13.00449
13.00449
13.00449
13.00449
13.00449
13.00449
TSDPT [45]
13.00495
13.00495
13.00495
13.00495
13.00495
13.00495
SSDPT [46]
13.00606
13.00606
13.00606
13.00606
13.00606
13.00606
Meshfree [37]
13.95287
13.95287
13.95287
13.95287
13.95287
13.95287
13.00516
13.00516
13.00516
13.00516
13.00516
13.00516
0.5
CLPT [37]
7.65398
8.25597
8.56223
8.78063
9.18254
9.61525
FSDPT [37]
7.33732
7.91320
8.20015
8.41034
8.78673
9.19517
TSDPT [45]
7.36437
7.94084
8.22470
8.43645
8.80997
9.21681
SSDPT [46]
7.36568
7.94195
8.22538
8.43712
8.81037
9.21670
Meshfree [37]
7.16207
7.71627
7.98956
8.22133
8.55172
8.94190
7.36361
7.94047
8.22048
8.43637
8.80980
9.21651
1
CLPT [37]
5.33248
6.02733
6.40391
6.68150
7.19663
7.78406
FSDPT [37]
5.14236
5.81379
6.17020
6.43892
6.92571
7.48365
TSDPT [45]
5.16713
5.84006
6.19394
6.46474
6.94944
7.50656
SSDPT [46]
5.16846
5.84119
6.19461
6.46539
6.94980
7.50629
Meshfree [37]
5.06137
5.71135
6.05467
6.31500
6.78405
7.31995
5.16489
5.83877
6.19209
6.48150
6.94876
7.50579
5
CLPT [37]
2.73080
3.10704
3.48418
3.65732
4.21238
4.85717
FSDPT [37]
2.63842
3.02252
3.38538
3.55958
4.09285
4.71475
TSDPT [45]
2.65821
3.04257
3.40351
3.57956
4.11209
4.73469
SSDPT [46]
2.66006
3.04406
3.40449
3.58063
4.11288
4.73488
Meshfree [37]
2.63652
3.00791
3.36225
3.53005
4.05070
4.64701
2.64204
3.02902
3.38599
3.57260
4.10499
4.73116
10
CLPT [37]
2.56985
2.80340
3.16427
3.25924
3.79238
4.38221
FSDPT [37]
2.46904
2.72626
3.07428
3.17521
3.68890
4.26040
TSDPT [45]
2.48727
2.74632
3.09190
3.19471
3.70752
4.27991
SSDPT [46]
2.48928
2.74844
3.13443
3.19456
3.14574
4.38175
Meshfree [37]
2.47216
2.72046
3.06067
3.15761
3.66166
4.20550
2.47035
2.72423
3.06339
3.18364
3.09531
4.27493Table 15
Biaxial buckling load of simply supported plate of C-type using the 5-order theory.
p
Theory
P
-
Source
1-0-1
2-1-2
2-1-1
1-1-1
2-2-1
1-2-1
0
CLPT [37]
6.86896
6.86896
6.86896
6.86896
6.86896
6.86896
FSDPT [37]
6.50224
6.50224
6.50224
6.50224
6.50224
6.50224
TSDPT [45]
6.50248
6.50248
6.50248
6.50248
6.50248
6.50248
SSDPT [46]
6.50303
6.50303
6.50303
6.50303
6.50303
6.50303
Meshfree [37]
6.47643
6.47643
6.47643
6.47643
6.47643
6.47643
6.50258
6.50258
6.50258
6.50258
6.50258
6.50258
0.5
CLPT [37]
3.82699
4.12798
4.28112
4.39032
4.59127
4.80762
FSDPT [37]
3.66866
3.95660
4.10007
4.20517
4.39336
4.59758
TSDPT [45]
3.68219
3.97042
4.11235
4.21823
4.40499
4.60841
SSDPT [46]
3.68284
3.97097
4.11269
4.21856
4.40519
4.60835
Meshfree [37]
3.58104
3.85813
3.99478
4.09639
4.27586
4.47095
3.68180
3.97024
4.11204
4.21819
4.40490
4.60826
1
CLPT [37]
2.66624
3.01366
3.20195
3.34075
3.59831
3.89203
FSDPT [37]
2.57118
2.90690
3.08510
3.21946
3.46286
3.74182
TSDPT [45]
2.58357
2.92003
3.09697
3.23237
3.47472
3.75328
SSDPT [46]
2.58423
2.92060
3.09731
3.23270
3.47490
3.75314
Meshfree [37]
2.53069
2.85568
3.02733
3.15750
3.39202
3.65998
2.58245
2.91939
3.09605
3.23212
3.47438
3.75290
5
CLPT [37]
1.36540
1.55352
1.74209
1.82866
2.10619
2.42859
FSDPT [37]
1.31921
1.51126
1.69269
1.77979
2.04642
2.35737
TSDPT [45]
1.32910
1.52129
1.70176
1.78978
2.05605
2.36734
SSDPT [46]
1.32910
1.52203
1.70224
1.79032
2.05644
2.36744
Meshfree [37]
1.31826
1.50395
1.68128
1.76502
2.02535
2.32351
1.32102
1.51451
1.69299
1.78630
2.05225
2.36558
10
CLPT [37]
1.28493
1.40170
1.58214
1.62962
1.89619
2.19111
FSDPT [37]
1.23452
1.36313
1.53714
1.58760
1.84445
2.13020
TSDPT [45]
1.24363
1.37316
1.54595
1.59736
1.85376
2.13995
SSDPT [46]
1.24475
1.37422
1.56721
1.59728
1.57287
2.19087
Meshfree [37]
1.23608
1.36023
1.53034
1.57880
1.83083
2.10275
1.23518
1.36211
1.53170
1.59182
1.84765
2.13746A good agreement between the present solution and references considered; especially [10, 29] is obtained. This allows us to conclude that the present higher-order plate theory is good for the modeling of simply supported sandwich FGM plates with OHOSNDPT as a good formulation. This study also lead us to conclude that the thickness stretching effect has a strong influence on the buckling analysis of sandwich FGM plates as gives higher fundamental buckling loads.The isotropic fully ceramic plate (first line on Tables14 and 15) has the higher fundamental buckling loads. As the core thickness to the total thickness of the plate ratio (
(
h
2
-
h
1
)
/
h
) increases, the buckling loads increase as well. Considering each column of both tables we may conclude that the critical buckling loads decrease as the power-law exponent p increases. By comparing Tables 14 and 15 we also conclude that the biaxial buckling load of simply supported sandwich square plate with FGM skins is half the uniaxial one for the same plate.


## 5.3.2. Isotropic FGM Plate Compared with Laminated Plate

In this subsection, the higher-order plate theories are used for the buckling analysis of simply supported functionally of two models of type D and E with side-to-thickness ratioa
/
h
=
10. All the material properties are taken as in Section 5.1.3.The convergence study with the higher-order theory was also conducted for each buckling load type considering the orderK of the plate theory from 1 to 7. The uniaxial case is presented in Table 16 for the Linear model with a
/
h
=
20,10,5. Further results are obtained by considering K equals 3, which seems acceptable by the convergence study.Table 16
Convergence study for the uniaxial buckling load of a simply supported square plate with linear model using the higher-order theory.
Scheme
K
=
1
K
=
2
K
=
3
K
=
4
K
=
5
K
=
6
K
=
7
a
/
h
=
20
30.6144
30.5958
30.5585
30.5584
30.5584
30.5584
30.5584
a
/
h
=
10
14.9028
14.8679
14.7979
14.7978
14.7978
14.7977
14.7977
a
/
h
=
5
6.7392
6.6849
6.5739
6.5739
6.5738
6.5738
6.5738In Table17, nondimensional uniaxial buckling load and biaxial buckling load of plate are computed for various a
/
h of three different models.Table 17
Bi and iniaxial buckling load of simply supported plate of type D using the higher-order theory.
Model
Reciprocal
Exponent
Linear
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
a
/
h
=
20
a
/
h
=
10
a
/
h
=
5
Biaxial
15.2792
7.3989
3.2869
10.9162
5.2696
2.3170
8.3944
4.0156
1.7136
Uniaxial
30.5585
14.7979
6.5739
21.8325
10.5393
4.6339
16.7888
8.0313
3.4273It is well known that as the plate thickness increases, the shear effect which leads to reduction of the critical buckling stress increases. It can be seen from Tables16 and 17 that although increasing the plate thickness decreases P
¯
x in all models, it can be concluded that for compressive biaxial and uniaxial loadings consideration of the prebuckling deformations in the stability equations increases the shear effect.Another interesting phenomenon appears between the static problem and the buckling situation. Whena
/
h
=
10, that is, the FGM is a thin plate, the rate of the buckling stresses between the reciprocal model and the linear model (the value is 1.8425) is almost the same as the proportion of the displacement at the center point (the more accurate value 1.8426 is than the number in Section 5.1.4) at the same K. When a
/
h
=
5, that is, the FGM is a moderately thick plate, the former ratio is 1.9181 and the latter one becomes 1.9198 (4.6322/2.4128). In summary, it not only illustrate that the present OHOSNDPT is high accuracy both the static problem and the buckling situation, but also reflect the stability for the thin plate and moderately thick plate.


## 6. Conclusions

Thus, in this paper, an attempt is made to study the analytical solutions of a special higher-order shear and normal deformable plate theory for static, free vibration and buckling analyses of thin and thick rectangular plate. Based on the special higher-order shear and normal deformable plate theory, the equations of motion are derived from Hamilton’s principle. The closed-form solution method and the algorithm for obtaining the numerical results are developed and employed. The numerical results on the displacement, stress distribution, the frequencies, and the critical buckling stress are presented and discussed. According to this discussion the following conclusions can be reached:(i)
The special higher-order shear and normal deformable plate theory uses various weighted orthogonal polynomial to solve arbitrary gradient FGM plate problems. A weighted orthogonal polynomial does correspond to only one material model. The own special polynomial has made the special higher-order shear and normal deformation theory, which is firstly named as the orthogonal higher-order shear and normal deformation theory (OHOSNDPT). The OHOSNDPT
(ii)
Hierarchical linear model with coupling OHOSNDPT is proposed. Compared with the hierarchical constant coefficient model coupled with OHOSNDPT and the common hierarchical constant coefficient model, the special hierarchical linear model is preferred to the two reference schemes. Furthermore, the two coupled schemes, whether linear model or constant coefficient model, are better than the common constant coefficient model.
(iii)
Finally, comparison studies are presented out for a large number of plates with different values of thickness ratio, multimaterial model, and various theory orderK. It is found that the distribution of the displacements and stress along the longitudinal coordinate-axis, the natural frequencies and buckling loads obtained by present theory match well with those obtained by the reference literatures.



---

*Source: 6879508-2017-01-03.xml*