Abstract
The isogeometric analysis with nonuniform rational B-spline (NURBS) based on the classical plate theory (CPT) is developed for free vibration analyses of functionally graded material (FGM) thin plates. The objective of this work is to provide an efficient and accurate numerical simulation approach for the nonhomogeneous thin plates and shells. Higher order basis functions can be easily obtained in IGA, thus the formulation of CPT based on the IGA can be simplified. For the FGM thin plates, material property gradient in the thickness direction is unsymmetrical about the midplane, so effects of midplane displacements cannot be ignored, whereas the CPT neglects midplane displacements. To eliminate the effects of midplane displacements without introducing new unknown variables, the physical neutral surface is introduced into the CPT. The approximation of the deflection field and the geometric description are performed by using the NURBS basis functions. Compared with the first-order shear deformation theory, the present method has lower memory consumption and higher efficiency. Several numerical results show that the present method yields highly accurate solutions.
1. Introduction
Although numerous plate theories have been developed since the late 19th century, the classical plate theory (CPT) [1] and the first-order shear deformation theory (FSDT) [2] are the most widely accepted and applied theories in engineering.
The CPT is the simplest theory with high computational efficiency that neglects transverse shear strains and midplane displacements. The effects of transverse shear strains and midplane displacements in homogeneous thin plates are small; hence, the CPT is widely used in homogeneous thin plates [3]. However, nonhomogeneous thin plates such as functionally graded material (FGM) plates have material property gradient in the thickness direction and this gradient is not symmetric about the midplane, so the results from the CPT may be incorrect. Our studies show that the natural frequencies of FGM thin plates based on the CPT and FSDT are close to each other for small or big gradient indices, while a significant difference exists for moderate gradient indices. As we all know, the transverse shear deformations have small effects in thin plates, so neglecting midplane displacements in the CPT causes differences in the FSDT solution for FGM thin plates. The material appears relatively homogeneous for small or big gradient indices; in other words, the material property gradient in the thickness direction is almost symmetric about the midplane, thus midplane displacements can be neglected, whereas the material is obviously nonhomogeneous for moderate gradient indices; this means that the material property gradient in the thickness direction is apparently unsymmetrical about the midplane, so neglecting midplane displacements in the CPT causes significant differences in the FSDT solution for FGM thin plates with moderate gradient indices. Abrate [4] indicated that stretching-bending couplings in constitutive equations of FGM plates do not exist if the proper reference surface is selected. Zhang and Zhou [5] recently proposed that the physical neutral surface is the proper reference surface and presented a theoretical analysis for FGM thin plates based on the physical neutral surface. To eliminate the effects of midplane displacements without introducing new unknown variables, we integrated the physical neutral surface into the CPT, that is, CPT_neu.
Regarding FGM plates of general geometries and loading conditions, analytical methods are usually infeasible, and numerical approaches are thus required. The finite element method is considered as one of the most versatile and powerful numerical methods. However, element-based polynomial approximation in the standard finite element method inevitably leads to discretization errors in complex geometry; the gap between the finite element analysis (FEA) and computer-aided design (CAD) consumes much time on the mesh generation of FEA. In 2005, Hughes et al. [6] developed a new numerical method called isogeometric analysis (IGA) to achieve a seamless integration between CAD and FEA. The principle of the IGA is that spline (such as NURBS and T-splines) basis functions are employed as shape functions for geometric description and field approximation. The unique advantages of the IGA are as follows [6]: (1) the same basis function space is adopted for geometric and mesh models, thus avoiding mesh discretization errors even on coarse mesh; (2) the natural division in the spline parametric domain is used instead of the mesh, thereby avoiding the traditional mesh generation procedure; (3) mesh refinement does not necessitate an external description of the geometry, thereby simplifying the adaptive mesh refinement; (4) the splines have higher order continuity, thus simplifying the construction of global higher order continuous conforming elements. In previous years, considerable work has been conducted in improving or applying the original IGA in some fields, including structural vibrations [7], plates and shells [8–10], fluid mechanics [11], fluid-structure interaction problems [12, 13], damage and fracture mechanics [14], and structural shape optimization [15].
The CPT needs higher order elements, so it is difficult to achieve it in tradition finite element method. In the IGA, higher order NURBS basis functions can be easily obtained, so the formulation of CPT can be simplified [1, 3]. The main objective of this research work is to extend CPT-based IGA to study free vibration of FGM thin plates. To serve this purpose, the CPT_neu is introduced into the NURBS-based IGA thus the method can analyze nonhomogeneous thin plates. Compared with the FSDT-based IGA, the present method has lower memory consumption and higher efficiency.
The paper is organized as follows. The NURBS-based IGA is described briefly in Section 2. In Section 3, the effects of midplane displacements on the natural frequency of FGM thin plates are investigated in detail, the physical neutral surface is then introduced into the CPT, and the CPT_neu-based eigenvalue equations of the IGA are developed. The numerical validation is presented in Section 4. The method is illustrated on several numerical examples in Section 5. Finally, conclusions are drawn in Section 6.
2. NURBS-Based Isogeometric Analysis
2.1. NURBS Basis Functions
NURBS basis functions are linear combinations of B-spline basis functions [16]. B-spline basis functions are used to parameterize curves and are generally defined at the knot vector
The knot vector
If w i (0 < w i ≤ 1) is the weight associated with the B-spline basis function N i, p (ξ), then the NURBS basis function R i, p (ξ), which is a weighted average of the B-spline basis functions, is defined as
The NURBS basis functions satisfy the partition of unity, that is,
Similarly, the two- or three-dimensional NURBS basis functions can be constructed by taking the tensor product of two or three one-dimensional B-spline basis functions:
where w
i, j
and w
i, j, k
denote the two- and three-dimensional geometry related weights, respectively; N
i, p
(ξ), N
j, q
(η), and N
k, r
(ζ) are the B-spline basis functions of order p in the ξ direction, order q in the η direction, and order r in the ζ direction, respectively; both N
j, q
(η) and N
k, r
(ζ) follow the recursive formula shown in (2) with knot vectors
By using the NURBS basis functions, a NURBS curve of order p can be constructed as
where
In a similar method, a two-dimensional NURBS surface and a three-dimensional NURBS solid can be constructed as
where
2.2. Isoparametric Discretization
Similar to the isoparametric finite element method, mapping from the parametric domain to the physical domain in the IGA is expressed as
where ξ is the parametric coordinate, that is, ξ in one dimension, (ξ, η) in two dimensions, and (ξ, η, ζ) in three dimensions;
Higher order basis functions can be easily obtained in the IGA, so only the deflection of plate w(
where w I is the deflection at control point I.
3. Fundamental Equations of FGM Thin Plate
3.1. FGM Plate
Consider a metal/ceramic FGM plate of thickness t. The metal is located at the bottom face of the plate, and the ceramic is located at the top face of the plate. The material is assumed to be isotropic, and the material properties are assumed to vary smoothly and continuously from the ceramic to the metal through the thickness. The xy-plane is the midplane of the plate with the z-axis in a positive upward from the midplane. Young's modulus E and mass density ρ are considered inhomogeneous material properties, whereas Poisson's ratio ν is assumed to be constant. The inhomogeneous material properties are expressed in power law form [18]:
where n is the volume fraction exponent or the gradient index; z is the thickness coordinate variable;
3.2. Effects of Midplane Displacements on FGM Thin Plate
The CPT based on the Kirchhoff hypothesis neglects transverse shear deformations and midplane displacements but can produce reasonably accurate results for homogeneous thin plates. To consider the effect of transverse shear deformation, shear deformation theories were proposed [2, 18]. For homogeneous thin plates, the accuracy of the CPT is almost equal to that of the shear deformation theories. Some studies show that the shear deformation theories can analyse FGM plates [19]. In the following section, we investigate whether the CPT can also effectively analyse FGM thin plates.
A simply supported aluminium/alumina (Al/Al2O3) square plate with length a = 1 m is considered. The plate is modelled with cubic NURBS elements. Length-to-thickness ratios of a/t = 20 and 100 are considered. A mesh with linear parameterization and 20 × 20 control points is used as shown in Figure 1. The material properties of Al and Al2O3 are v m = 0.3, E m = 70 GPa, and ρ m = 2707 kg/m3 and v c = 0.3, E c = 380 GPa, and ρ m = 3800 kg/m3, respectively. Figure 2 presents a comparison of the first normalized free vibration frequency of the plate obtained by the NURBS-based IGA based on the CPT and FSDT. The solutions based on the CPT and the FSDT are close to each other for small or big gradient indices. A significant difference exists between the two theories for moderate gradient indices. The length-to-thickness ratio has insignificant effects on the difference between the CPT and the FSDT for the FGM thin plate, and this means that transverse shear deformations can be neglected in thin plate. The material appears relatively homogeneous for small or big gradient indices, whereas the material is obviously nonhomogeneous for moderate gradient indices. Differences exist between the CPT and the FSDT because the FSDT considers the effects of transverse shear deformations and midplane displacements. The previous analysis shows that neglecting the midplane displacements caused by inhomogeneous materials in the CPT may cause significant differences in the FSDT solution for moderate gradient indices.

A full square plate with 20 × 20 control points and 17 × 17 elements: (a) control mesh, (b) physical mesh.

Comparison of the first normalized frequency of a simply supported Al/Al2O3 thin plate based on the CPT and the FSDT: (a) a/t = 20,
3.3. FGM Thin Plate Equations Based on the Physical Neutral Surface
Considering a FGM rectangular plate under pure bending with the theory of elasticity, the coordinate system and the material parameter variation are the same as those in Section 3.1. One surface (z = z0) whose normal strains and stresses are zero can be found, and the surface is called as physical neutral surface [5]. The definition of physical neutral surface is expressed as [5]
From (12), the physical neutral surface and geometric middle surface are the same for homogeneous isotropic plates.
Using the physical neutral surface instead of geometry neutral surface in the CPT, the effects of midplane displacements can be avoided. To eliminate the effects of the midplane displacements in the FGM thin plate without introducing new unknown variables, the physical neutral surface concept is integrated into the CPT, that is, CPT_neu. By using the physical neutral surface and the CPT, the displacement fields can be assumed as the following:
The strains are defines as
According to the generalized Hooke's law, the stresses can be obtained as
where
The pseudostrains and pseudostresses of the plate are denoted as follows:
where
For the free vibration analyses, the Galerkin weak form of the elastodynamic undamped equilibrium equation can be expressed as follows [20]:
Substituting (10), (17) into (19) introduces the dynamic discrete system of equations:
where
The element contribution to
and the element contribution to
with
The nth normal mode φ n can be obtained by separating the variables:
where ω n is the nth natural frequency.
Substituting (24) into (20) produces the eigenvalue equations of the FGM plates:
4. Numerical Validation
To validate the accuracy of the proposed method, the Al/Al2O3 thin plate considered in the previous section is again employed. All conditions are the same as those in the previous section. Figure 3 shows the comparison of the first normalized natural frequencies of the plate obtained by the proposed method and the NURBS-based IGA based on the FSDT. For the length-to-thickness ratios of a/t = 20 and 100, the results obtained by the NURBS-based IGA based on the CPT_neu match well with those from the NURBS-based IGA based on the FSDT. It is concluded that the effects of midplane displacements vanish in the CPT_neu.

Comparison of the first normalized frequency of a simply supported Al/Al2O3 thin plate based on the CPT_neu and the FSDT: (a) a/t = 20,
5. Numerical Results
In this section, several examples with reference solutions are provided to assess the robustness of the proposed method. All examples reported later are calculated by using the cubic order NURBS basis functions with linear parameterization [6]. A 4 × 4 Gauss quadrature scheme is used in each NURBS element. The same Al/Al2O3 material considered in the previous section is employed.
5.1. A Simply Supported Al/Al2O3 Square Thin Plate with Length a = 1 m
The length-to-thickness ratios of a/t = 20 and 100 are considered. The natural frequencies of the plates are investigated, and the results are compared with the NURBS-based IGA based on the FSDT or analytical solutions [17, 21]. Figure 4 shows the comparison of the first normalized natural frequencies from the proposed method and the reference solutions with different control points. The following features are evident: (1) results based on the CPT_neu and FSDT match well with the analytical solutions [17, 21]; (2) as the CPT, natural frequencies obtained by the CPT_neu are larger than the analytical solutions; (3) the accuracy of the CPT_neu is higher than that of the FSDT when the number of control points is small, whereas the accuracy of the CPT_neu is almost the same as that of the FSDT with increasing number of control points.

Comparison of the first normalized frequency of a simply supported Al/Al2O3 thin plate based on the CPT_neu, FSDT, and analytical solutions: (a) a/t = 20,
The computational time versus the number of control points is plotted in Figure 5. The results indicate that the CPT_neu is much more efficient than the FSDT because five degrees of freedom exist in one point for the FSDT, whereas only one degree of freedom exists in one point for the CPT_neu.

Comparison of the computational time between CPT_neu and FSDT.
5.2. An Al/Al2O3 Rectangular Thin Plate with a Length a and Width b under Different Boundary Conditions
The effects of the gradient indices, length-to-width ratios (a/b), and boundary conditions of the FGM thin plates on the natural frequencies are investigated to demonstrate the validity of the present method. For this purpose, the normalized natural frequencies of two length-to-width ratio plates with four gradient indices and three boundary conditions are examined; the results are then compared with the analytical solutions [17] and the NURBS-based IGA based on the CPT and the FSDT. The three boundary conditions are as follows: (1) fully simply supported (SSSS); (2) left and right boundaries are simply supported, whereas bottom and upper boundaries are free (SFSF); (3) left and right boundaries are simply supported, whereas bottom and upper boundaries are clamped (SCSC). The plate length is a = 1 m, and a mesh with 20 × 20 control points is used.
Tables 1 and 2 present the first five modes of the normalized frequencies for a/b = 1 and a/b = 2, respectively. The exact solutions [17] and the results obtained by the numerical approaches based on the CPT and the FSDT are also presented for comparison. The frequencies computed by the present method match well with the exact solutions and the solutions by the NURBS-based IGA based on the FSDT. The results from the present method are better than those obtained from the NURBS-based IGA based on the CPT.
Normalized natural frequencies
Normalized natural frequencies
6. Conclusions
Higher order basis functions can be easily obtained in the IGA, thus the formulation of CPT based on the IGA can be simplified. For nonhomogeneous thin plate such as FGM, the midplane displacements are neglected in the CPT, so the CPT is not suitable for analyzing obviously nonhomogeneous thin plate such as FGM thin plates with moderate gradient indices. To eliminate the effects of midplane displacements without introducing new unknown variables, the physical neutral surface is integrated into the CPT, that is, CPT_neu. The NURBS-based isogeometric analysis with the CPT_neu is developed to solve natural frequencies of FGM thin plates. Numerical simulations indicate that the accuracy of the CPT_neu is almost the same as that of the FSDT. The CPT_neu has one degree of freedom in one point, whereas the FSDT has five degrees of freedom in one point. Hence, the result of the CPT_neu is much more efficient than that of the FSDT.
Footnotes
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant no. 51179063) and the Fundamental Research Funds for the Central Universities (Grant no. 2011B03014) and Jiangsu Province Graduate Students Research and Innovation Plan (Grant no. CXZZ13_023).
