Static and Free Vibration Analysis of Laminated Composite Plates Using Isogeometric Approach Based on the Third Order Shear Deformation Theory

Abstract

Isogeometric analysis (IGA) based on nonuniform rational B-splines (NURBS) is applied for static and free vibration analysis of laminated composite plates by using the third order shear deformation theory (TSDT). TSDT requires C¹-continuity of generalized displacements and NURBS basis functions are well-suited for this requirement. Due to the noninterpolatory nature of NURBS basis functions, a penalty method is applied to enforce the essential boundary conditions. The validity and accuracy of the present method are demonstrated through a series of numerical examples of isotropic and laminated composite plates with different shapes, boundary conditions, fiber orientations, lay-up numbers, and so forth. The obtained numerical results are compared with either the analytical solutions or other available numerical methods.

1. Introduction

The requirements for higher stiffness and strength-to-weight ratio, better wear resistance, and greater health characteristics over conventional materials have resulted in the wide application of composite laminates in the area of aerospace, automotive, and underwater structures. Therefore, static and dynamic analysis of both thin and thick laminates becomes an important research topic. Based on different assumptions for displacement fields, several theories for analysis of composite laminates have been developed. The classical plate theory (CPT) relying on the Kirchhoff-Love assumptions [1, 2] can only be applied if plates are thin, because transverse shear deformations are not considered. The first order shear deformation theory (FSDT) which accounts for transverse shear effects [3, 4] can be applied for both moderately thick and thin plates. However, in FSDT, the transverse shear strain is assumed to be constant along the thickness direction, and, thus a shear correction factor (SCF) has to be considered. The SCF depends on many factors, such as material coefficients, stacking scheme, plate geometry and boundary conditions; therefore the evaluation of the SCF remains a subject of research [5]. In addition, with constant shear strain assumption, transverse shear stresses obtained may not be zero on the plate surface, which contradicts the traction free boundary conditions.

To overcome the limitations of FSDT, a more accurate high order shear deformation theory, termed the third order shear deformation theory (TSDT) [6, 7], has been developed to formulate the mechanical behavior of composite laminates. TSDT introduces the cubic variation of the displacement such that more accurate interlaminar stress distributions can be yielded. Therefore, the SCF does not need to be calculated, and the conditions of zero shear stress at both the top and the bottom plate surfaces can be automatically satisfied. However, TSDT requires C¹-continuity of generalized displacements, which is not an easy task to achieve using the C⁰-continuity based finite element method (FEM). To overcome this difficulty, several methods including some mesh-free approaches have been used, for example, the moving least-squares (MLS) method [8], the radial point interpolation method (RPIM) [9], the radial basis function collocation method [10], the finite point formulation with the multiquadric radial basis function method [11], and so forth. In this paper, we demonstrate that such a C¹-TSDT formulation can be conveniently achieved by using NURBS-based isogeometric analysis.

Isogeometric analysis (IGA) was first proposed by Hughes and coworkers [12, 13]. The basic idea is that basis functions, usually NURBS basis functions, which accurately represent the geometry are directly used as the interpolation functions of the unknown field variables. Thus, it offers a possibility to close the existing gaps between computer aided design (CAD) and finite element analysis (FEA). Exact geometries are fixed at the coarsest level of discretization and mesh refinement process can be vastly simplified by eliminating the necessity to communicate with the CAD description of the geometry. There are three methods of mesh refinement in IGA, that is, h-refinement of knot insertion, p-refinement of order elevation, and the unique k-refinement of increasing the continuity of basis functions. Therefore, higher order NURBS basis functions with increased interelement continuity can be directly obtained.

Since the pioneer work of Bazilevs and Hughes described above, IGA has been widely applied to various engineering problems including fluid mechanics [14], fluid-structure interaction problems [15, 16], structural shape optimization [17, 18], shell analysis [19–21], linear and nonlinear elasticity and plasticity behavior [22, 23], structural vibration [24], and so forth. In the field of plate analysis employing NURBS-based IGA, Shojaee et al. [25] conducted the analysis of thin composite laminates based on the CPT; Thai et al. [26] performed the analysis of composite laminates using the FSDT; Tran et al. [27] first studied the behavior of functionally graded material plates based on the TSDT. However, numerical examples of these previous researches usually confined themselves into plates with simple shapes and the effects of three plate theories (i.e., CPT, FSDT, and TSDT) on numerical results have not been fully discussed yet.

This paper further develops the NURBS-based isogeometric approach on static and free vibration analysis of composite laminates by using the TSDT. Different from the previous researches, primary focuses are placed on exploring the advantages of the present method in the analysis of plates with complicated shapes. The effects of three plate theories on numerical results are discussed as well.

The remainder of the paper is organized as follows: first, B-Splines, NURBS, and NURBS-based IGA are introduced in Section 2; according to TSDT, an isogeometric formulation for laminated plates is derived in Section 3; Section 4 presents several numerical examples for both static and free vibration analyses of composite laminates with different physical parameters and boundary conditions. Numerical results are compared with analytical solutions and the available data obtained by other numerical methods. Finally, some concluding remarks are given in Section 5.

2. NURBS and NURBS-Based IGA

2.1. Knot Vectors and B-Splines Basis Functions

A knot vector, Ξ, is a set of nondecreasing coordinates in the parametric space

$\mathrm{\Xi }=\left\{{\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\dots ,{\xi }_{n+p+\mathrm{1}}\right\},$

(1)

where ξ _i ∈ [0, 1] is the ith knot, i is the knot index, i = 1, 2,…, n + p + 1, p is the order of the B-spline, and n is the number of basis functions. A knot vector is said to be open if its first and last knots are repeated (p + 1) times [28]. Open knot vectors are standard in the CAD literature. Unless otherwise stated, open knot vectors are always used in this paper.

With the given knot vector Ξ, B-spline basis functions can be constructed recursively by the Cox-de Boor recursion formula [13]; it starts from p = 0 with

$N_{i,\mathrm{0}}\left(\xi \right)=\{\begin{array}{cc}\mathrm{0}\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}{\xi }_i\le \xi <{\xi }_{i+\mathrm{1}}\hfill \\ \mathrm{1}\hfill & \text{otherwise},\hfill \end{array}$

(2)

and for p = 1, 2, 3,… they are defined by

$N_{i,p}\left(\xi \right)=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-\mathrm{1}}\left(\xi \right)+\frac{{\xi }_{i+p+\mathrm{1}}-\xi }{{\xi }_{i+p+\mathrm{1}}-{\xi }_{i+\mathrm{1}}}{N}_{i+\mathrm{1},p-\mathrm{1}}\left(\xi \right).$

(3)

2.2. Nonuniform Rational B-Splines (NURBS)

Non-uniform rational B-splines are called NURBS for short. NURBS basis function is given by

$R_i^p\left(\xi \right)=\frac{N_{i,p}\left(\xi \right){\omega }_i}{W\left(\xi \right)}=\frac{N_{i,p}\left(\xi \right){\omega }_i}{{\sum }_{\widehat{i}=\mathrm{1}}^n{N}_{\widehat{i},p}\left(\xi \right){\omega }_{\widehat{i}}},$

(4)

which means that every B-spline basis function N _i,p (ξ) is set to have a weight ω _i [13].

A NURBS curve [13] can be constructed by using (4) in conjunction with the control points B _i ,

$\mathbf{C}\left(\xi \right)=\sum _{i=\mathrm{1}}^n{R}_i^p\left(\xi \right){\mathbf{B}}_i.$

(5)

Given two knot vectors Ξ = {ξ₁, ξ₂,…, ξ_{n + p + 1}}, H = {η₁, η₂,…, η_{m + q + 1}} and control points B _i,j , a two-dimensional NURBS surface [13] can be defined as

$\mathbf{S}\left(\xi ,\eta \right)=\sum _{i=\mathrm{1}}^n\hspace{0.17em}\sum _{j=\mathrm{1}}^m{R}_{i,j}^{p,q}\left(\xi ,\eta \right){\mathbf{B}}_{i,j},$

(6)

where R _i,j ^p,q (ξ, η) is the tensor product of NURBS basis functions of two parametric dimensions ξ and η,

$R_{i,j}^{p,q}\left(\xi ,\eta \right)=\frac{N_{i,p}\left(\xi \right)M_{j,q}\left(\eta \right){\omega }_{i,j}}{{\sum }_{\widehat{i}=\mathrm{1}}^n{\sum }_{\widehat{j}=\mathrm{1}}^m{N}_{\widehat{i},p}\left(\xi \right)M_{\widehat{j},q}\left(\eta \right){\omega }_{\widehat{i},\widehat{j}}}.$

(7)

An example of quadratic NURBS surface with its control points represented by red dots is illustrated in Figure 1.

Figure 1:

Quadratic NURBS surface with its control points.

2.3. NURBS-Based IGA

IGA utilizes the concept of isoparametric, which refers to the use of the same basis functions for both geometry and unknown field discretization. The displacements of the control points are the degrees of freedom of the structure. NURBS patch defines a subdomain and knot spans subdivide the subdomain into elements. Mesh refinement strategies are developed from a combination of knot insertion and order elevation, resulting in the h, p or k-refinements as mentioned above. The support of each basis function consists of a small number of elements, depending on the polynomial degree and multiplicity of the knot. Numerical integration is done by Gauss quadrature formulation at the element level.

3. An Isogeometric Formulation for Composite Laminates Based on the TSDT

3.1. The Third Order Shear Deformation Theory

Consider a laminated plate of thickness h shown in Figure 2. The x-y plane is taken to be the midplane of the laminate before deformation, and the kth layer is located between z(k) and z(k + 1) in the thickness direction.

Figure 2:

A typical laminated plate and its coordinate system.

The displacement field of a plate based on the TSDT is as follows [7],

$\begin{array}{c}u=u_{\mathrm{0}}+z{\phi }_x-\alpha z^{\mathrm{3}}\left({\phi }_x+\frac{\partial w}{\partial x}\right)\\ v=v_{\mathrm{0}}+z{\phi }_y-\alpha z^{\mathrm{3}}\left({\phi }_y+\frac{\partial w}{\partial y}\right).\\ w=w_{\mathrm{0}},\end{array}$

(8)

where α = 4/(3 h²); (u₀, v₀, w₀) are the displacements of a point on the reference plane in the (x, y, z) directions, respectively; (φ _x , φ _y ) are the rotations with respect to y and x axis as shown in Figure 2.

While TSDT is applied, five degrees of freedom, u₀, v₀, w₀, φ _x , and φ _y are considered. Equation (8) can also be expressed in a matrix form as

$\left\{\begin{array}{c}u\\ v\\ w\end{array}\right\}\hspace{0.17em}=\left[\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& -\alpha z^{\mathrm{3}}\frac{\partial }{\partial x}& z-z^{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& -\alpha z^{\mathrm{3}}\frac{\partial }{\partial y}& \mathrm{0}& z-\alpha z^{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right]\left\{\begin{array}{c}{u}_{\mathrm{0}}\\ v_{\mathrm{0}}\\ w_{\mathrm{0}}\\ {\phi }_x\\ {\phi }_y\end{array}\right\},\hspace{0.5em}\text{or}\hspace{1em}\mathbf{u}=\mathbf{H}{\mathbf{u}}_{\mathrm{0}}.$

(9)

It can be noted that by setting α = 0, the displacement field based on the first order shear deformation theory (FSDT) can be obtained. Furthermore, by defining α = 0, φ _x = − ∂w/∂x, and φ _y = − ∂w/∂y, the displacement field based on the classical plate theory (CPT) can be obtained [29].

The pseudostrains with respect to the displacement field are

$\left[\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\\ {\epsilon }_{xz}\\ {\epsilon }_{yz}\end{array}\right]=\left[\begin{array}{c}\frac{\partial u}{\partial x}\\ \frac{\partial v}{\partial y}\\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\\ \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\\ \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\end{array}\right]=\left[\begin{array}{ccccc}\frac{\partial }{\partial x}& \mathrm{0}& -\alpha z^{\mathrm{3}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}& \left(z-\alpha z^{\mathrm{3}}\right)\frac{\partial }{\partial x}& \mathrm{0}\\ \mathrm{0}& \frac{\partial }{\partial y}& -\alpha z^{\mathrm{3}}\frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}& \mathrm{0}& \left(z-\alpha z^{\mathrm{3}}\right)\frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}& -\mathrm{2}\alpha z^{\mathrm{3}}\frac{{\partial }^{\mathrm{2}}}{\partial x\partial y}& \left(z-\alpha z^{\mathrm{3}}\right)\frac{\partial }{\partial y}& \left(z-\alpha z^{\mathrm{3}}\right)\frac{\partial }{\partial x}\\ \mathrm{0}& \mathrm{0}& \left(\mathrm{1}-\beta z^{\mathrm{2}}\right)\frac{\partial }{\partial x}& \left(\mathrm{1}-\beta z^{\mathrm{2}}\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \left(\mathrm{1}-\beta z^{\mathrm{2}}\right)\frac{\partial }{\partial y}& \mathrm{0}& \left(\mathrm{1}-\beta z^{\mathrm{2}}\right)\end{array}\right]\left\{\begin{array}{c}{u}_{\mathrm{0}}\\ v_{\mathrm{0}}\\ w_{\mathrm{0}}\\ {\phi }_x\\ {\phi }_y\end{array}\right\},$

(10)

where β = 3α. Equation (10) can also be stated in a more compact form as

${\mathbf{\epsilon }}_p=\mathbf{L}{\mathbf{u}}_{\mathrm{0}},$

(11)

where L is the differential operator matrix.

The strain-stress relations [30] for an orthotropic lamina in the plate coordinate system (x, y, z) are

$\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\\ {\sigma }_{xz}\\ {\sigma }_{yz}\end{array}\right\}\hspace{0.5em}=\left[\begin{array}{ccccc}{\stackrel{-}{Q}}_{\mathrm{11}}& {\stackrel{-}{Q}}_{\mathrm{12}}& {\stackrel{-}{Q}}_{\mathrm{16}}& \mathrm{0}& \mathrm{0}\\ {\stackrel{-}{Q}}_{\mathrm{21}}& {\stackrel{-}{Q}}_{\mathrm{22}}& {\stackrel{-}{Q}}_{\mathrm{26}}& \mathrm{0}& \mathrm{0}\\ {\stackrel{-}{Q}}_{\mathrm{61}}& {\stackrel{-}{Q}}_{\mathrm{62}}& {\stackrel{-}{Q}}_{\mathrm{66}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\stackrel{-}{Q}}_{\mathrm{44}}& {\stackrel{-}{Q}}_{\mathrm{45}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\stackrel{-}{Q}}_{\mathrm{54}}& {\stackrel{-}{Q}}_{\mathrm{55}}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\\ {\epsilon }_{xz}\\ {\epsilon }_{yz}\end{array}\right\}\hspace{1em}\text{or}\hspace{0.17em}\hspace{0.17em}{\mathit{\sigma }}_p=\mathbf{D}{\mathbf{\epsilon }}_p,$

(12)

where

$\begin{array}{c}{\stackrel{-}{Q}}_{\mathrm{11}}=Q_{\mathrm{11}}{\cos }^{\mathrm{4}}\theta +\mathrm{2}\left(Q_{\mathrm{12}}+\mathrm{2}{Q}_{\mathrm{66}}\right){\sin }^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta +Q_{\mathrm{22}}{\sin }^{\mathrm{4}}\theta ,\\ {\stackrel{-}{Q}}_{\mathrm{12}}={\stackrel{-}{Q}}_{\mathrm{21}}=\left(Q_{\mathrm{11}}+Q_{\mathrm{22}}-\mathrm{4}{Q}_{\mathrm{66}}\right){\sin }^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta +Q_{\mathrm{12}}\left({\sin }^{\mathrm{4}}\theta +{\cos }^{\mathrm{4}}\theta \right),\\ {\stackrel{-}{Q}}_{\mathrm{22}}=Q_{\mathrm{11}}{\sin }^{\mathrm{4}}\theta +\mathrm{2}\left(Q_{\mathrm{12}}+\mathrm{2}{Q}_{\mathrm{66}}\right){\sin }^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta +Q_{\mathrm{22}}{\cos }^{\mathrm{4}}\theta ,\\ {\stackrel{-}{Q}}_{\mathrm{16}}={\stackrel{-}{Q}}_{\mathrm{61}}=\left(Q_{\mathrm{11}}-Q_{\mathrm{12}}-\mathrm{2}{Q}_{\mathrm{66}}\right)\sin \theta {\cos }^{\mathrm{3}}\theta +\left(Q_{\mathrm{12}}-Q_{\mathrm{22}}+\mathrm{2}{Q}_{\mathrm{66}}\right)\text{si}{\text{n}}^{\mathrm{3}}\theta \cos \theta ,\\ {\stackrel{-}{Q}}_{\mathrm{26}}={\stackrel{-}{Q}}_{\mathrm{62}}=\left(Q_{\mathrm{11}}-Q_{\mathrm{12}}-\mathrm{2}{Q}_{\mathrm{66}}\right)\text{si}{\text{n}}^{\mathrm{3}}\theta \cos \theta +\left(Q_{\mathrm{12}}-Q_{\mathrm{22}}+\mathrm{2}{Q}_{\mathrm{66}}\right)\sin \theta {\cos }^{\mathrm{3}}\theta ,\\ {\stackrel{-}{Q}}_{\mathrm{66}}=\left(Q_{\mathrm{11}}+Q_{\mathrm{22}}-\mathrm{2}{Q}_{\mathrm{12}}-\mathrm{2}{Q}_{\mathrm{66}}\right)\text{si}{\text{n}}^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta +Q_{\mathrm{66}}\left({\sin }^{\mathrm{4}}\theta +{\cos }^{\mathrm{4}}\theta \right),\\ {\stackrel{-}{Q}}_{\mathrm{44}}=Q_{\mathrm{44}}{\cos }^{\mathrm{2}}\theta +Q_{\mathrm{55}}{\sin }^{\mathrm{2}}\theta ,\\ {\stackrel{-}{Q}}_{\mathrm{45}}=\left(Q_{\mathrm{55}}-Q_{\mathrm{44}}\right)\cos \theta \sin \theta ,\\ {\stackrel{-}{Q}}_{\mathrm{55}}=Q_{\mathrm{55}}{\cos }^{\mathrm{2}}\theta +Q_{\mathrm{44}}{\sin }^{\mathrm{2}}\theta .\end{array}$

(13)

All Q _ij in (13) are defined in the material coordinate system of the lamina and θ is the angle of fiber orientation. Q _ij are given by

$\begin{array}{c}{Q}_{\mathrm{11}}=\frac{E_{\mathrm{1}}}{\mathrm{1}-{\nu }_{\mathrm{12}}{\nu }_{\mathrm{21}}},\hspace{1em}\hspace{1em}{Q}_{\mathrm{12}}=\frac{{\nu }_{\mathrm{12}}{E}_{\mathrm{2}}}{\mathrm{1}-{\nu }_{\mathrm{12}}{\nu }_{\mathrm{21}}},\hspace{1em}\hspace{1em}\\ Q_{\mathrm{22}}=\frac{E_{\mathrm{2}}}{\mathrm{1}-{\nu }_{\mathrm{12}}{\nu }_{\mathrm{21}}},\hspace{1em}\hspace{1em}{Q}_{\mathrm{66}}=G_{\mathrm{12}},\\ Q_{\mathrm{44}}=G_{\mathrm{13}},\hspace{1em}\hspace{1em}{Q}_{\mathrm{55}}=G_{\mathrm{23}},\end{array}$

(14)

where the subscript 1 refers to the fiber direction and subscript 2 refers to the direction transverse to the fiber [30].E₁ and E₂ are the Young modulus in the 1 and 2 directions, respectively; G₁₂, G₂₃, and G₁₃ are the shear modulus in the 1–2, 2–3, and 3–1 planes, respectively; ν _ij are Poisson ratios.

3.2. Penalty Method and Variational Form of System Equations

The boundary conditions for a clamped boundary of a plate can be denoted as follows:

$\begin{array}{c}{u}_{\mathrm{0}}=\mathrm{0},\hspace{1em}\hspace{1em}{v}_{\mathrm{0}}=\mathrm{0},\hspace{1em}\hspace{1em}{w}_{\mathrm{0}}=\mathrm{0},\\ {\phi }_{\mathbf{n}}=\mathrm{0},\hspace{1em}\hspace{1em}{\phi }_{\mathbf{s}}=\mathrm{0},\hspace{1em}\hspace{1em}\frac{\partial w}{\partial \mathbf{n}}=\mathrm{0},\end{array}$

(15)

where n = (n _x , n _y ) is the unit vector normal to the essential boundary and s is the unit vector perpendicular to n . ${\phi }_{\mathit{n}}$ and ${\phi }_{\mathit{s}}$ can be decomposed by the following equations:

${\phi }_{\mathbf{n}}=n_x{\phi }_x+n_y{\phi }_y,\hspace{1em}\hspace{1em}{\phi }_{\mathbf{s}}=n_x{\phi }_y-n_y{\phi }_x.$

(16)

Therefore, (15) can be rearranged as

$\left[\begin{array}{c}{u}_{\mathrm{0}}\\ v_{\mathrm{0}}\\ w_{\mathrm{0}}\\ {\phi }_{\mathit{n}}\\ {\phi }_{\mathit{s}}\\ \frac{\partial w}{\partial \mathit{n}}\end{array}\right]=\left[\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& n_x& n_y\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -n_y& n_x\\ \mathrm{0}& \mathrm{0}& \frac{\partial }{\partial \mathbf{n}}& \mathrm{0}& \mathrm{0}\end{array}\right]\left\{\begin{array}{c}{u}_{\mathrm{0}}\\ v_{\mathrm{0}}\\ w_{\mathrm{0}}\\ {\phi }_x\\ {\phi }_y\end{array}\right\}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right],\hspace{1em}\text{or}\hspace{0.17em}\hspace{0.17em}\tilde{\mathbf{u}}={\mathbf{L}}_{\mathbf{b}}{\mathbf{u}}_{\mathrm{0}}=\mathbf{0},$

(17)

where ${\mathit{L}}_{\mathit{b}}$ is the differential operator matrix.

The boundary conditions for a simply supported boundary are

$u_{\mathrm{0}}=\mathrm{0},\hspace{1em}\hspace{1em}{v}_{\mathrm{0}}=\mathrm{0},\hspace{1em}\hspace{1em}{w}_{\mathrm{0}}=\mathrm{0},\hspace{1em}\hspace{1em}{\phi }_{\mathbf{s}}=\mathrm{0}.$

(18)

Therefore, the differential operator matrix ${\mathit{L}}_{\mathit{b}}$ for a simply supported boundary is

${\mathbf{L}}_{\mathbf{b}}=\left[\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& 1& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& 1& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -n_y& n_x\end{array}\right].$

(19)

Since NURBS basis functions do not satisfy the Kronecker-Delta property, special treatments of essential boundary conditions have to be employed using the ideas in meshless method. Available techniques include Lagrange multiplier method [31], penalty method [32], Nitsche's method [33], and others [34]. In this study, penalty method is applied to enforce the essential boundary conditions in both static deflection analysis and free vibration analysis.

Considering the penalty method, the variational form of the equilibrium equation can be expanded as

${\int }_S\delta {\mathbf{\epsilon }}_p^T:{\mathit{\sigma }}_p\hspace{0.17em}\text{d}V+{\int }_V\delta {\mathbf{u}}^T·\rho \ddot{\mathbf{u}}\hspace{0.17em}\text{d}V-{\int }_V\delta {\mathbf{u}}^T·\mathbf{b}\hspace{0.17em}\text{d}V-{\int }_{{\mathrm{\Gamma }}_{\sigma }}\delta {\mathbf{u}}^T·\stackrel{-}{\mathbf{t}}\hspace{0.17em}\text{d}\mathrm{\Gamma }-\delta {\int }_{{\mathrm{\Gamma }}_u}{\left(\tilde{\mathbf{u}}-\stackrel{-}{\mathbf{u}}\right)}^T·\mathit{\alpha }·\left(\tilde{\mathbf{u}}-\stackrel{-}{\mathbf{u}}\right)\text{d}\mathrm{\Gamma }=\mathbf{0},$

(20)

where V and Γ are, respectively, the plate volume and plate surface, (σ _p ,ε _p ) are the pseudostress and pseudostrain, ρ is the density, b is the body force, $\stackrel{-}{\mathbf{t}}$ is the surface force, and α is a diagonal matrix of the penalty coefficient. In this study, the penalty coefficient is set to 1 × 10¹⁰ [8].

Substituting (9), (11), (12) and (17) into (20), the variational form of the equilibrium equation can be rewritten as

${\int }_S\delta {\mathbf{u}}_{\mathrm{0}}^T{\mathbf{L}}^T\mathbf{D}\mathbf{L}{\mathbf{u}}_{\mathrm{0}}\hspace{0.17em}\text{d}V+{\int }_V\rho \delta {\mathbf{u}}_{\mathrm{0}}^T{\mathbf{H}}^T\mathbf{H}{\ddot{\mathbf{u}}}_{\mathrm{0}}\hspace{0.17em}\text{d}V-{\int }_V\delta {\mathbf{u}}_{\mathrm{0}}^T{\mathbf{H}}^T\mathbf{b}\hspace{0.17em}\text{d}V-{\int }_{{\mathrm{\Gamma }}_{\sigma }}\delta {\mathbf{u}}_{\mathrm{0}}^T{\mathbf{H}}^T\stackrel{-}{\mathbf{t}}\hspace{0.17em}\text{d}\mathrm{\Gamma }-\delta {\int }_{{\mathrm{\Gamma }}_u}{\left({\mathbf{L}}_{\mathbf{b}}{\mathbf{u}}_{\mathrm{0}}-\stackrel{-}{\mathbf{u}}\right)}^T\mathit{\alpha }\left({\mathbf{L}}_{\mathbf{b}}{\mathbf{u}}_{\mathrm{0}}-\stackrel{-}{\mathbf{u}}\right)\text{d}\mathrm{\Gamma }=\mathbf{0}.$

(21)

3.3. Discrete System Equations

Using NURBS basis functions, the displacement field can be interpolated as

${\mathbf{u}}_{\mathrm{0}}=\sum _{I=\mathrm{1}}^n{R}_I{\mathbf{u}}_{\mathrm{0}}^I,$

(22)

where R _I are the NURBS basis functions and n is total number of them; ${\mathbf{u}}_{\mathrm{0}}^I={[u_{\mathrm{0}}^I\hspace{0.17em}\hspace{0.17em}{v}_{\mathrm{0}}^I\hspace{0.17em}\hspace{0.17em}{w}_{\mathrm{0}}^I\hspace{0.17em}\hspace{0.17em}{\phi }_x^I\hspace{0.17em}\hspace{0.17em}{\phi }_y^I]}^T$ is the vector of nodal degrees of freedom associated with each control point.

Substituting (22) into the variational form (21), the final discrete system equations can be obtained as follows:

$\mathbf{M}{\ddot{\mathbf{u}}}_{\mathrm{0}}+\left(\mathbf{K}+\tilde{\mathbf{K}}\right){\mathbf{u}}_{\mathrm{0}}=\mathbf{F},$

(23)

where stiffness matrix K , $\tilde{\mathbf{K}}$ and mass matrix M are given by

$\begin{array}{c}{\mathbf{K}}_{IJ}={\int }_S\left(\sum _{i=\mathrm{1}}^{N_L}{\int }_{z_i}{\mathbf{B}}_{Ii}^T{\mathbf{D}}_i{\mathbf{B}}_{Ji}\text{d}z\right)\text{d}S,\\ {\tilde{\mathbf{K}}}_{IJ}={\int }_S\left(\sum _{i=\mathrm{1}}^{N_L}{\int }_{z_i}{\tilde{\mathbf{B}}}_{Ii}^T\mathit{\alpha }{\tilde{\mathbf{B}}}_{Ji}\text{d}z\right)\text{d}S,\\ {\mathbf{M}}_{IJ}={\int }_S\left(\sum _{i=\mathrm{1}}^{N_L}{\int }_{z_i}{\rho }_i{\mathbf{N}}_{Ii}^T{\mathbf{N}}_{Ji}\text{d}z\right)\text{d}S,\end{array}$

(24)

where N _L denotes the number of layers of the laminated plate. During computation, stiffness and mass matrices should be firstly calculated at each layer and then assembled into system matrix together.

B _Ii has the form B _Ii = L R _I , where L is the differential operator matrix shown in (10). Similarly, ${\tilde{\mathbf{B}}}_{Ii}={\mathit{L}}_{\mathit{b}}{R}_I$ [29].

N _Ii is expressed as [9]

${\mathbf{N}}_{Ii}=\left[\begin{array}{ccccc}{R}_I& \mathrm{0}& -\alpha z^{\mathrm{3}}{R}_{I,x}& \left(z-\alpha z^{\mathrm{3}}\right)R_I& \mathrm{0}\\ \mathrm{0}& R_I& -\alpha z^{\mathrm{3}}{R}_{I,y}& \mathrm{0}& \left(z-\alpha z^{\mathrm{3}}\right)R_I\\ \mathrm{0}& \mathrm{0}& R_I& \mathrm{0}& \mathrm{0}\end{array}\right].$

(25)

In the static analysis of this paper, only transverse pressure p applied on the plate surface is considered; the global force vector F is simplified as

$\mathbf{F}={\int }_{S}p{R}_I\text{d}S.$

(26)

For static analysis, the formulation of a laminated plate can be expressed as

$\left(\mathbf{K}+\tilde{\mathbf{K}}\right)\mathbf{U}=\mathbf{F},$

(27)

and for free vibration analysis, the formulation can be written as

$\left[\left(\mathbf{K}+\tilde{\mathbf{K}}\right)-{\omega }^{\mathrm{2}}\mathbf{M}\right]\mathbf{U}=\mathrm{0},$

(28)

where ω is the natural frequency.

4. Numerical Examples

In this section, static and free vibration behaviors of both isotropic and laminated composite plates are studied based on the TSDT. Simple square and rectangular plates to more complicated shape plates with various aspect ratios and different boundary conditions are considered. Unless otherwise stated, layers in composite laminates are assumed to have identical thickness and material properties, except that fiber orientations are different. Numerical results are obtained using full (p + 1) × (q + 1) Gauss integration points, where p and q are the orders of the NURBS basis functions in two parametric directions ξ and η, respectively.

4.1. Static Deflection Analysis

4.1.1. Isotropic Elliptical Plate under a Uniform Transverse Pressure

To demonstrate the validity and accuracy of IGA in modeling common engineering shapes, a fully clamped elliptical plate under a uniform transverse pressure is chosen to be analyzed. The lengths of the major and minor axes are assumed to be a = 2 m and b = 1 m; the aspect ratio a/h is taken as 100. A uniform load P = − 1 N is applied. The material of the plate is isotropic with elastic modulus E = 10⁴ MPa and Poisson ratio v = 0.3. Control points and physical meshes of 10 × 10 and 15 × 15 uniform meshed elliptical plates modeled by quadratic NURBS basis functions are shown in Figure 3, respectively.

Figure 3:

Control points and physical meshes of (a) 10 × 10 and (b) 15 × 15 uniform meshed elliptical plates.

A series of uniform meshes with quadratic, cubic, and quartic NURBS elements are used in this numerical study. The obtained center displacements are normalized by analytical solutions [35] and the results are plotted against the number of elements per side in Figure 4. It can be seen that the present results converge to the exact solutions. The quartic NURBS elements give sufficiently good results even when the mesh is 10 × 10; however, quadratic NURBS elements suffer from shear locking and cannot get reasonable results on coarse meshes.

Figure 4:

Normalized center displacements versus the number of elements per side of a fully clamped elliptical plate.

4.1.2. Four Layers Cross-Ply Square Laminated Plate (0°/90°/90°/0°) under a Sinusoidal Distributed Transverse Load

A four-layer cross-ply square laminated plate (0°/90°/90°/0°) under a sinusoidal distributed transverse load q is now considered as shown in Figure 5. The simply supported boundary condition is applied with load q defined as

$q\left(x,y\right)=-q_{\mathrm{0}}\cos \frac{\pi }{a}\left(x-\frac{a}{\mathrm{2}}\right)\sin \frac{\pi }{a}y.$

(29)

Figure 5:

A square laminated plate under a sinusoidal distributed load.

The material properties are given by E₁ = 25E₂, G₁₂ = G₁₃ = 0.5E₂, G₂₃ = 0.2E₂, and ν₁₂ = 0.25. By employing the cubic NURBS basis functions, different physical mesh densities are used to study the convergence of the present method. In this numerical study, three densities with 2 × 2, 5 × 5, and 10 × 10 uniform meshes are considered. Physical meshes and control points of the three models are shown in Figure 6. It is clearly evident that geometries can remain precise during h-refinement process. Four plate aspect ratios ranging from 4 to 100, representing thick to relatively thin plates, are employed for simulation. The results of the center deflections are presented using the following nondimensional form for all aspect ratio cases:

$\stackrel{-}{w}=w\left(\frac{a}{\mathrm{2}},\frac{a}{\mathrm{2}}\right)\frac{h^{\mathrm{3}}{E}_{\mathrm{2}}}{a^{\mathrm{4}}{q}_{\mathrm{0}}}\times \mathrm{1}{\mathrm{0}}^{\mathrm{2}}.$

(30)

Figure 6:

Physical meshes and control points of a square plate with (a) 2 × 2 uniform mesh, (b) 5 × 5 uniform mesh, and (c) 10 × 10 uniform mesh.

Nondimensional center deflections obtained by the three given densities of physical meshes are presented in Table 1. In order to compare TSDT with other plate theories, computations based on CPT and FSDT are also conducted on exactly the same mesh employing isogeometric analysis. In addition, other available solutions obtained by mesh-free method of Dai et al. [8], finite point formulation with the multiquadric radial basis function of Ferreira et al. [11], FEM of Reddy [30], and 3D elasticity theory (ELS) of Pagano and Hatfield [36] are also given in Table 1.

Table 1:

Nondimensional center deflections of a square laminated plate under a sinusoidal distributed transverse load.

Plate theory	a/h
	4	10	20	100
Elasticity
Pagano and Hatfield [36]	1.954	0.743	0.517	0.4347
TSDT
2 × 2 mesh	1.9052	0.7198	0.5079	0.4118
5 × 5 mesh	1.8928	0.7144	0.5058	0.4336
10 × 10 mesh	1.8937	0.7148	0.5060	0.4343
Dai et al. [8]	1.8930	0.7147	0.5060	0.4342
Ferreira et al. [11]	1.8864	0.7153	0.5070	0.4365
Reddy [30]	1.894	0.715	0.506	0.438
FSDT
2 × 2 mesh	1.7218	0.6658	0.4911	0.4018
5 × 5 mesh	1.7088	0.6624	0.4910	0.4328
10 × 10 mesh	1.7096	0.6627	0.4912	0.4337
Dai et al. [8]	1.7310	0.6708	0.4950	0.4340
Reddy [30]	1.710	0.663	0.491	0.434
CPT
2 × 2 mesh	0.4348
5 × 5 mesh	0.4310
10 × 10 mesh	0.4313
Dai et al. [8]	0.4302
Reddy [30]	0.431

It can be seen from Table 1 that isogeometric analysis gives solutions that are very agreeable to those of FEM and mesh-free methods. Numerical results can be improved by h-refinement especially when the laminates are thin. The differences of the results between FSDT and ELS get larger when the aspect ratio decreases, for example, when a/h = 100, the result of FSDT is 0.16% smaller than that of ELS; however, when a/h = 4, the relative difference reaches 12.51%, whereas for TSDT, the difference remains very small as 3.09%. This means that FSDT is not suitable for the analysis of thick laminates. It is also noted that the CPT results do not change with aspect ratios, simply because the CPT does not consider thickness effect of laminates. Therefore, just as expected, TSDT is more accurate than FSDT and CPT for not only thin laminates but also thick ones.

In order to investigate the suitability of the present method for the simulation of stresses in composite laminates, stress analysis is conducted in this numerical example. The normalized stresses are given by

$\begin{array}{c}{\stackrel{-}{\sigma }}_{xx}=\frac{h^{\mathrm{2}}}{q_{\mathrm{0}}{a}^{\mathrm{2}}}{\sigma }_{xx}\left(\frac{a}{\mathrm{2}},\frac{a}{\mathrm{2}},\frac{h}{\mathrm{2}}\right),\hspace{1em}\hspace{1em}{\stackrel{-}{\sigma }}_{yy}=\frac{h^{\mathrm{2}}}{q_{\mathrm{0}}{a}^{\mathrm{2}}}{\sigma }_{yy}\left(\frac{a}{\mathrm{2}},\frac{a}{\mathrm{2}},\frac{h}{\mathrm{4}}\right);\\ {\stackrel{-}{\tau }}_{xz}=\frac{h}{q_{\mathrm{0}}b}{\tau }_{xz}\left(\mathrm{0},\frac{a}{\mathrm{2}},\mathrm{0}\right),\hspace{1em}\hspace{1em}{\stackrel{-}{\tau }}_{xy}=\frac{h^{\mathrm{2}}}{q_{\mathrm{0}}{b}^{\mathrm{2}}}{\tau }_{xy}\left(\mathrm{0,0},\frac{h}{\mathrm{2}}\right).\end{array}$

(31)

By employing cubic NURBS basis functions and 10 × 10 meshes, the normalized stresses of the present method are shown in Table 2. The results of several other methods including the elasticity solution given by Pagano [37], the finite strip method (FSM) based on FSDT by Akhras et al. [38], the finite element method (FEM) based on TSDT proposed by Reddy [7], and the multiquadric radial basis functions (RBFs) method based on TSDT reported by Ferreira et al. [11] are also listed in Table 2. It can be seen that the present solutions agree very well with other published results, especially the numerical results of Reddy [7].

Table 2:

Normalized stresses of a square laminated plate under a sinusoidal distributed transverse load.

a/h	Method	${\stackrel{-}{\sigma }}_{xx}$	${\stackrel{-}{\sigma }}_{yy}$	${\stackrel{-}{\tau }}_{zx}$	${\stackrel{-}{\tau }}_{xy}$
4	Pagano [37]	0.720	0.666	0.270	0.0467
	Akhras et al. [38]	0.4059	0.5765	0.1398	0.0308
	Reddy [7]	0.6651	0.6322	0.2064	0.0440
	Ferreira et al. [11]	0.6659	0.6313	0.1352	0.0433
	Present	0.6624	0.6319	0.2064	0.0440
10	Pagano [37]	0.559	0.403	0.301	0.0276
	Akhras et al. [38]	0.4989	0.3615	0.1667	0.0241
	Reddy [7]	0.5456	0.3888	0.2640	0.0268
	Ferreira et al. [11]	0.5466	0.4383	0.3347	0.0267
	Present	0.5446	0.3887	0.2640	0.0268
20	Pagano [37]	0.543	0.309	0.328	0.0230
	Akhras et al. [38]	0.5273	0.2957	0.1749	0.0221
	Reddy [7]	0.5393	0.3043	0.2825	0.0228
	Ferreira et al. [11]	0.5405	0.3648	0.3818	0.0228
	Present	0.5386	0.3042	0.2825	0.0228
100	Pagano [37]	0.539	0.271	0.339	0.0214
	Akhras et al. [38]	0.5382	0.2705	0.1780	0.0213
	Reddy [7]	0.5387	0.2708	0.2897	0.0213
	Ferreira et al. [11]	0.5413	0.3359	0.4106	0.0215
	Present	0.5381	0.2708	0.2898	0.0213

4.1.3. Three Layers Square Sandwich Plate under a Uniform Transverse Pressure

A simply supported square sandwich plate [39] subjected to a uniform transverse load q is considered in this numerical example.

The material properties of the sandwich core are expressed in the stiffness matrix, ${\stackrel{-}{Q}}_{\text{core}}$ as

$\begin{array}{c}\begin{array}{c}{\stackrel{-}{Q}}_{\text{core}}=\left[\begin{array}{ccccc}\mathrm{0.999781}& \mathrm{0.231192}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0.231192}& \mathrm{0.524886}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0.262931}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0.159914}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0.266810}\end{array}\right].\end{array}\end{array}$

(32)

Face material properties are related to the core material properties as

${\stackrel{-}{Q}}_{\text{face}}=R{\stackrel{-}{Q}}_{\text{core}}.$

(33)

The aspect ratio a/h is taken as 10. A 10 × 10 uniform mesh with cubic NURBS basis functions is used in this numerical study. Normalized displacements and stresses of the square sandwich plates are defined as

$\begin{array}{c}\stackrel{-}{\omega }=\omega \frac{\mathrm{0.999781}}{hq},\hspace{1em}\hspace{1em}{\stackrel{-}{\sigma }}_x^{\mathrm{1}}=\frac{{\sigma }_x^{\left(\mathrm{3}\right)}\left(a/\mathrm{2},a/\mathrm{2},h/\mathrm{2}\right)}{q},\\ {\stackrel{-}{\sigma }}_x^{\mathrm{2}}=\frac{{\sigma }_x^{\left(\mathrm{3}\right)}\left(a/\mathrm{2},a/\mathrm{2},\mathrm{2}h/\mathrm{5}\right)}{q},\hspace{1em}\hspace{1em}{\stackrel{-}{\sigma }}_x^{\mathrm{3}}=\frac{{\sigma }_x^{\left(\mathrm{2}\right)}\left(a/\mathrm{2},a/\mathrm{2},\mathrm{2}h/\mathrm{5}\right)}{q},\\ {\stackrel{-}{\sigma }}_y^{\mathrm{1}}=\frac{{\sigma }_y^{\left(\mathrm{3}\right)}\left(a/\mathrm{2},a/\mathrm{2},h/\mathrm{2}\right)}{q},\hspace{1em}\hspace{1em}{\stackrel{-}{\sigma }}_y^{\mathrm{2}}=\frac{{\sigma }_y^{\left(\mathrm{3}\right)}\left(a/\mathrm{2},a/\mathrm{2},\mathrm{2}h/\mathrm{5}\right)}{q},\\ {\stackrel{-}{\sigma }}_y^{\mathrm{3}}=\frac{{\sigma }_y^{\left(\mathrm{2}\right)}\left(a/\mathrm{2},a/\mathrm{2},\mathrm{2}h/\mathrm{5}\right)}{q},\hspace{1em}\hspace{1em}{\stackrel{-}{\tau }}_{xz}^{\mathrm{1}}=\frac{{\tau }_{xz}^{\left(\mathrm{2}\right)}\left(\mathrm{0},a/\mathrm{2},\mathrm{0}\right)}{q},\\ {\stackrel{-}{\tau }}_{xz}^{\mathrm{2}}=\frac{{\tau }_{xz}^{\left(\mathrm{2}\right)}\left(\mathrm{0},a/\mathrm{2},-\mathrm{2}h/\mathrm{5}\right)}{q},\end{array}$

(34)

where superscript 1, 2, and 3 during the parentheses are layer numbers.

As shown in Table 3, the obtained results with different R values are compared with available solutions by the FEM method of Pandya and Kant [40] and multiquadric radial basis function method (RBF) reported by Ferreira et al. [11]. These two published results are also based on TSDT. In addition, exact solutions of Srinivas [39] are also provided as a reference.

Table 3:

Normalized displacements and stresses of the square sandwich plates.

R	Method	$\stackrel{-}{\omega }$	${\stackrel{-}{\sigma }}_x^1$	${\stackrel{-}{\sigma }}_x^2$	${\stackrel{-}{\sigma }}_x^3$	${\stackrel{-}{\sigma }}_y^1$	${\stackrel{-}{\sigma }}_y^2$	${\stackrel{-}{\sigma }}_y^3$	${\stackrel{-}{\tau }}_{xz}^1$	${\stackrel{-}{\tau }}_{xz}^2$
5	Srinivas [39]	258.97	60.353	46.623	9.340	38.491	30.097	6.161	4.364	3.268
	Pandya and Kant [40]	256.13	62.380	46.910	9.382	38.930	30.330	6.065	3.089	2.566
	Ferreira et al. [11]	257.110	60.366	47.003	9.401	38.456	30.242	6.048	4.548	2.391
	Present	256.969	60.320	46.972	9.394	38.424	30.219	6.044	4.776	2.942
10	Srinivas [39]	159.38	65.332	48.857	4.903	43.566	33.413	3.500	4.096	3.515
	Pandya and Kant [40]	152.33	64.650	51.310	5.131	42.830	33.970	3.397	3.147	2.587
	Ferreira et al. [11]	154.658	65.381	49.973	4.997	43.240	33.637	3.364	3.528	2.398
	Present	154.499	65.310	49.923	4.992	43.190	33.599	3.360	4.312	2.656
15	Srinivas [39]	121.72	66.787	48.299	3.238	46.424	34.955	2.494	3.964	3.577
	Pandya and Kant [40]	110.43	66.620	51.970	3.465	44.920	35.410	2.361	3.035	2.691
	Ferreira et al. [11]	114.644	66.920	50.323	3.355	45.623	35.170	2.345	3.021	2.275
	Present	114.497	66.838	50.264	3.351	45.561	35.123	2.342	3.961	2.440

It can be seen that for all three R values considered, normal stresses σ _x and σ _y obtained by the present method are in excellent agreement with the exact solutions and the results of the other two numerical methods. As for shear stresses τ _xz , the present method can produce reasonable results at location z = 0, which is in layer 2; however, at location z = − 2 h/5, which is the interlamination interface between layers 1 and 2, all these three methods cannot produce accurate results. This is probably due to the discontinuity and the sudden change of material properties at the interlamination interface.

The stress distributions along the thickness of the sandwich plate with R = 15 are plotted in Figure 7 for normal stress component σ _x and shear stress component τ _xz , respectively. As shown in Figure 7(a), due to the plate being thin, both TSDT and FSDT can predict the normal stress reasonably except that TSDT can get more accurate results than FSDT near the plate surface. Figure 7(b) gives the shear stress distributions and it can be seen that TSDT predicts zero shear stress at both the bottom and the top surfaces, while the results of FSDT remain constant within each layer and are not zero at the plate surfaces, which contradict the real situation.

Figure 7:

The distribution of (a) normal stress σ _x and (b) shear stress τ _xz along thickness.

4.2. Natural Frequency Analysis

4.2.1. Clamped L-Shaped Isotropic Plate

A clamped L-shaped isotropic plate with thickness to length ratio h/a = 0.2 and Poisson ratio v = 0.3 is considered in this numerical example [43]. The geometry and NURBS meshes are provided in Figure 8. The natural frequencies are normalized by

$\stackrel{-}{\omega }=\omega a\sqrt{\frac{\rho }{G}},$

(35)

where G = E/2(1 + v). In Table 4, nondimensional frequency parameters of the present method are compared with the results obtained by the finite element method using reduced integration Q4 plate elements based on the FSDT by Ferreira [41], the global collocation RBF technique based on TSDT also reported by Ferreira et al. [42], the local RBF finite difference technique based on TSDT provided by Roque et al. [43], and the IGA method based on FSDT given by Thai et al. [26]. It can be seen that with the increase of the order of NURBS elements and the improvement of mesh quality, numerical results of the present method converge to the results of the other published methods. The first four mode shapes of the L-shaped plate obtained by the present method are shown in Figure 9.

Table 4:

Nondimensional frequency parameters of a clamped isotropic L-shaped plate.

Method	Mesh	${\stackrel{-}{\omega }}_{\mathrm{1}}$	${\stackrel{-}{\omega }}_{\mathrm{2}}$	${\stackrel{-}{\omega }}_{\mathrm{3}}$	${\stackrel{-}{\omega }}_{\mathrm{4}}$
Quadratic	10 × 5	2.6161	3.1529	3.6436	4.8434
	20 × 10	2.0279	2.5892	2.9655	3.8563
	30 × 15	1.9356	2.4906	2.8576	3.7259
Cubic	10 × 5	2.0608	2.6169	3.0040	3.9227
	20 × 10	1.9044	2.4560	2.8181	3.6801
	30 × 15	1.8678	2.4171	2.7752	3.6341
Quartic	10 × 5	1.9465	2.5000	2.8630	3.7293
	20 × 10	1.8663	2.4154	2.7727	3.6320
	30 × 15	1.8454	2.3937	2.7474	3.6069
FEM [41]		1.8557	2.3952	2.7637	3.6210
Global [42]		1.8468	2.3266	2.6806	3.4987
Local [43]		1.8832	2.3450	2.7698	3.5714
IGA-FSDT [26]		1.8147	2.3741	2.7570	3.6091

Figure 8:

(a) Geometry, (b) 10 × 5, (c) 20 × 10, and (d) 30 × 15 NURBS meshes of a clamped L-shaped isotropic plate.

Figure 9:

The first four mode shapes of a clamped L-shaped plate.

4.2.2. Isotropic Square Plate with a Complicated Cutout

In order to demonstrate the applicability of isogeometric analysis to exactly represent complicated shapes, natural frequencies are computed for an isotropic square plate with a complicated cutout shown in Figure 10(a). As shown in Figure 10(b), the plate is divided into 8 NURBS patches [31]. The material and geometrical properties of the plate are given by elastic modulus E = 2 × 10¹¹ MPa, Poisson ratio v = 0.3, mass density ρ = 8000 kg/m³, and plate thickness h = 0.05 m.

Figure 10:

(a) Geometric parameters (unit: m) and (b) discretization of the patches of an isotropic square plate with a complicated cutout.

With only one element in each patch, the exact geometry can be represented with only eight elements using 28 control points as shown in Figure 11(a). By dividing each patch with 5 × 4 uniform mesh using h-refinement, followed by increasing the order of NURBS basis functions using k-refinement, the refined mesh can be obtained for the numerical study. Figure 11(b) shows the physical mesh and control points of the complicated shape plate modeled by cubic NURBS basis functions.

Figure 11:

(a) The coarsest mesh and (b) refined mesh of an isotopic square plate with a hole of complicated shape. Blue lines mark the element boundaries, black lines represent the patch boundaries, and red circles are control points.

Natural frequencies are normalized by $\stackrel{-}{\omega }={({\omega }^2{a}^4\rho h/D)}^{1/4}$ , where a is the plate edge length and D is the bending stiffness of the plate defined as D = Eh³/(12(1 − v²)). The first ten normalized natural frequencies obtained by quadratic, cubic, and quartic NURBS basis functions are listed in Tables 5 and 6 for simply supported and clamped boundary conditions, respectively. Numerical solutions of the present method are compared with available results obtained by isogeometric analysis (IGA) based on the classical plate theory (CPT) [31], moving Kriging interpolation-based mesh-free method (MKI) [44], element free Galerkin method (EFG) [45], and radial point interpolation method (RPIM) [46], and also solutions obtained by commercial FEM software ABAQUS. It can be observed that results of the present method are in good agreement with solutions of RPIM and ABAQUS for both boundary conditions while the results of MKI and EFG methods are generally a little higher. It can be noted that the results of the present method have slight differences with the IGA-CPT solutions [31] this can be explained by using different plate theories. Just as with the case in static analysis, CPT does not consider thickness effect of plates. Convergence of numerical results is observed during k-refinement process. The first six mode shapes of a simply supported plate obtained by the present method are shown in Figure 12.

Table 5:

Normalized natural frequencies of a simply supported isotropic square plate with a complicated cutout.

Quadratic	Cubic	Quartic	IGA-CPT [31]	MKI [44]	EFG [45]	RPIM [46]	ABAQUS
4.959	4.913	4.911	5.193	5.390	5.453	4.919	4.910
6.447	6.388	6.382	6.579	7.502	8.069	6.398	6.376
6.875	6.760	6.749	6.597	8.347	9.554	6.775	6.741
8.762	8.563	8.551	7.819	10.636	10.099	8.613	8.547
9.192	8.978	8.958	8.812	11.048	11.328	9.016	8.942
11.102	10.673	10.641	9.42	12.894	12.765	10.738	10.629
11.488	10.922	10.880	10.742	13.710	13.685	10.930	10.873
13.080	11.673	11.600	10.776	14.062	14.305	11.601	11.584
13.523	12.839	12.796	11.919	16.694	15.721	12.903	12.790
14.064	13.210	13.142	13.200	17.364	17.079	13.283	13.122

Table 6:

Normalized natural frequencies of a clamped isotropic square plate with a complicated cutout.

Quadratic	Cubic	Quartic	IGA-CPT [31]	EFG [45]	RPIM [46]	ABAQUS
8.954	7.726	7.558	7.621	7.548	7.410	7.428
12.137	10.352	10.056	9.810	10.764	9.726	9.809
12.338	10.438	10.106	9.948	11.113	9.764	9.830
13.219	11.540	11.191	11.135	11.328	10.896	10.929
13.405	11.639	11.361	11.216	12.862	11.114	11.130
14.576	12.854	12.529	12.482	13.300	12.353	12.302
15.404	13.566	13.110	12.872	14.168	12.781	12.812
16.958	14.761	13.763	13.650	15.369	13.368	13.416
17.342	15.056	14.645	14.676	16.205	14.485	14.408
19.561	15.299	14.913	14.738	17.137	14.766	14.660

Figure 12:

The first six mode shapes of a simply supported isotropic square plate with a complicated cutout.

4.2.3. Four Layers Cross-Ply Square Laminated Plate (0°/90°/90°/0°) with a Complicated Cutout

A four-layer cross-ply square laminated plate (0°/90°/90°/0°) with a complicated cutout [47] is analyzed as the last numerical example of this section. The geometry of the plate is shown in Figure 13(a). As shown in Figure 13(b), the plate is divided into 9 NURBS patches. The material properties for all layers of the laminate are identical: E₁ = 40E₂, G₂₃ = 0.5E₂, G₁₃ = G₁₂ = 0.6E₂, and ν₁₂ = 0.25. The original coarsest mesh with only nine elements and 34 control points is shown in Figure 14(a). Using the same refinement algorithm as that in Section 4.2.1, Figure 14(b) gives the physical mesh and control points of the refined mesh modeled by quadratic NURBS basis functions. A thick laminated plate with an aspect ratio of 5 and an extremely thin laminated plate with an aspect ratio of 1000 are analyzed, respectively. Only the simply supported boundary condition is considered for simplicity. The first ten normalized natural frequencies are listed in Tables 7 and 8 for the two different aspect ratios, respectively. FEM results obtained by ABAQUS are also provided for comparison.

Table 7:

Normalized natural frequencies of a square laminate with a complicated cutout of aspect ratio 5.

Quadratic	Cubic	Quartic	ABAQUS
5.710	5.582	5.571	5.638
6.788	6.687	6.678	6.988
7.851	7.747	7.738	7.200
11.799	11.638	11.623	12.124
13.786	13.606	13.591	13.686
16.324	16.122	16.109	15.965
17.180	16.899	16.876	16.840
19.808	19.402	19.363	19.699
20.493	20.069	20.035	19.805
23.097	22.673	22.627	22.042

Table 8:

Normalized natural frequencies of a square laminate with a complicated cutout of aspect ratio 1000.

Quadratic	Cubic	Quartic	ABAQUS
11.077	10.035	9.978	8.600
17.272	14.986	14.924	13.627
25.999	23.030	22.790	22.149
41.368	36.598	36.155	36.167
51.074	38.175	37.978	37.461
73.370	49.699	48.801	48.336
74.722	66.729	65.701	65.936
101.441	77.333	74.903	76.404
105.636	84.106	82.465	82.415
137.553	98.360	94.272	94.992

Figure 13:

(a) Geometric parameters (Unit: m) and (b) discretization of the patches of a square laminated plate with a complicated cutout.

Figure 14:

(a) The coarsest mesh and (b) refined mesh of a square laminated plate with a complicated cutout.

From Tables 7 and 8 it can be seen that isogeometric analysis based on TSDT performs well when the plates are thick, even quadratic NURBS basis functions can get reasonable results; however, when the plates become extremely thin, low order NURBS basis functions suffer from shear locking. When the order increases, numerical results can be improved. The results of quartic NURBS basis functions are in good agreement with those of ABAQUS. The first six mode shapes of a simply supported square laminate with aspect ratio 1000 obtained by the present method are provided in Figure 15.

Figure 15:

The first six mode shapes of a simply supported square laminate with aspect ratio 1000.

5. Conclusion

A NURBS-based isogeometric approach has been applied for static and free vibration analysis of laminated composite plates using the third order shear deformation theory. Several isotropic and laminated composite plates with different shapes and boundary conditions are examined for various aspect ratios from thick to extremely thin scenarios. Comparisons are made with the numerical results obtained by other existing methods; isogeometric analysis based on TSDT better matches the actual situation and yields very good results, especially when the plates are thick. One of the advantages of IGA is increasing the order of basis functions without any difficulty. Therefore, high order NURBS basis functions are developed to improve the numerical results especially when the plates are extremely thin. The other advantage of IGA lies in the exact descriptions of complicated geometries. Numerical results of this study demonstrate that the present method is accurate and robust in performing free vibration analyses of plates with complicated shapes.

In the present approach, the phenomenon of shear locking has been noticed and it can be alleviated by mesh refinement. However, the computational cost will be increased due to an excessive increment of control points and the high order of NURBS basis functions. Locking free numerical methods [48–50] and parallel computations [51] are already successfully implemented in isogeometric analysis. Our further studies will focus on their applications in the analysis of beams, plates and shells.

Conflict of Interests

All the authors listed have reviewed the final version of the paper and approved to submit it to this journal. The authors do not have any conflict of interests regarding the content of this paper.