Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

53, 2, pp. 453-466, Warsaw 2015
DOI: 10.15632/jtam-pl.53.2.453

FREE VIBRATION AND BUCKLING ANALYSIS OF HIGHER ORDER
LAMINATED COMPOSITE PLATES USING THE ISOGEOMETRIC

APPROACH

Ognjen Peković, Slobodan Stupar, Aleksandar Simonović,

Jelena Svorcan, Srd-an Trivković

University of Belgrade, Faculty of Mechanical Engineering, Serbia

e-mail: opekovic@mas.bg.ac.rs

This research paper presents a higher order isogeometric laminated composite plate finite
element formulation. The isogeometric formulation is based onNURBS (non-uniform ratio-
nal B-splines) basis functions of varying degree. Plate kinematics is based on the third order
shear deformation theory (TSDT) of Reddy in order to avoid shear locking. Free vibration
and the buckling response of laminated composite plates are obtained and efficiency of the
method is considered. Numerical results with different element order are presented and the
obtained results are compared to analytical and conventional numerical results as well as
existing isogeometric plate finite elements.

Keywords: isogeometric, laminated composite plate, third order shear deformation theory

1. Introduction

Laminated composites arewidely used in aerospace,marine andwind turbine industries. Recen-
tly, there has appeared a great number of general industrial products made of composites. The
reasons for this are high strength-to-weight ratio, high stiffness, good dimensional stability after
manufacturing and high impact, fatigue and corrosion resistance of composites. In addition to
this, composite laminates possess ability to follow complex mould shapes and to be specifically
tailored through optimization of ply numbers and fibre orientations through the structure so
that they canmeet specific needs while minimizing weight (Jones, 1999).

Laminates in general have thickness much smaller than their planar dimensions so one can
use various plate theories instead of general 3D elasticity equations for their analysis (Reddy,
2004). Laminatedplate theories are classified into three groups: 1) equivalent single layer theories
(ESL), 2) layerwise plate theories and 3) individual-layer plate theories (Nosier et al., 1993). The
equivalent single layer theories are most widely used because they provide sufficiently accurate
description of the global laminate response at low computational cost. AmongESL, the classical
plate theory (CPT) is the simplest, but gives accurate results only for very thin plates since it
is unable to capture transverse shear effects. The first order shear deformation theories (FSDT)
give constant through thickness transverse shear strains resulting in constant transverse shear
stresses through thickness, which is contradictory to the elasticity solution. In order tomake up
for this, onemust use shear correction factors that are hard to determine. Higher order theories
introduce additional unknown variables but are capable of modelling realistic transverse shear
strains. Among them, the third order shear deformation theory of Reddy (1984) is the most
popular among engineers. It uses quadratic variation of transverse shear strains with vanishing
transverse shear stresses on the top and bottom surface, thus eliminating the need for the use
of shear correction factors.



454 O. Peković et al.

For arbitrary shapes and boundary conditions, the governing plate equations cannot be
solved analytically. Among different numerical techniques that seek approximate solutions, the
finite element method became a standard tool for treatment of stress analysis problems. In
FEM, the unknown field variables are approximated by linear combination of interpolation
(trial or shape) functions. In the standard FEM formulation, interpolation functions are locally
definedpolynomials inside theelement andzero everywhere outside theconsidered element.Most
existing finite elements and all commercial codes use Lagrangian (C0 interelement continuity)
and Hermitian (C1 interelement continuity) basis functions.
A great need exists in industry for integration of themanufacturing process from conceptual

phaseanddesign (bymeans of computer-aideddesign (CAD)) throughanalysis (byusingcompu-
ter aided engineering (CAE) tools) to manufacturing (done on CNCmachines trough computer
aided manufacturing (CAM)). CAD and CAM industries rely on the use of NURBS geometry
(Piegl, 1997; Rogers, 2001) for shape representation, thus CAD/CAM integration is relatively
straightforward. Although specialized CAD/CAM/CAE systems exist for the last 20-25 years
(PTC Creo, CATIAV5...) communication between CAD and CAE is not straightforward, and
it is necessary to build a new finite element model in order to run the necessary analysis. This
task takes up to 80% of total analysis time and is therefore one of the major bottlenecks in
CAD/CAE/CAM integration (Cottrell et al., 2009).
In order to overcome those difficulties, a new technique formally named isogeometric analysis

is proposed (Hughes et al., 2005). It allows the execution of analysis on geometrical CADmodel.
InsteadofLagrangeorHermitpolynomialbasis functions, the isogeometric finiteelementmethod
relies on NURBS basis functions, same as almost every CAD or CAM package. NURBS offers
generalmathematical representation of bothanalytical geometric objects and freeformgeometry.
Recently, several research papers have appeared that used the isogeometric approach for plate
and shell analysis (Kiendl et al., 2009; Benson et al., 2011; Echter et al., 2013) and composite
plate and shell analysis (Shojaee et al., 2012; Thai et al., 2012, 2013; Casanova and Gallego,
2013).
This paper presents free vibration and buckling analysis of TSDT composite plates in isoge-

ometric framework. Isogeometric formulations of global stiffness, mass and geometrical stiffness
matrix for quadratic, cubic and quartic elements are presented. All global matrices are formu-
lated in full compliance with the TSDT of Reddy. Results are compared to other available data
to demonstrate the accuracy of the proposedmethod.

2. NURBS curves, basis functions and surfaces

Non-uniform rational B-spline (NURBS) can represent arbitrary freeform shapeswith analytical
exactness that is needed in computer graphics (CG),CADandCAMapplications. After decades
of technology improvement, NURBS provides users with great control over the object shape in
an intuitive way with low memory consumption making them the most widespread technology
for shape representation.
NURBS are generalizations of nonrational Bezier and nonrational B-splines curves and sur-

faces. Bezier curves are parametric polynomial curves defined as

C(u)=
n∑

i=0

Bi,n(u)Pi 0¬ u ¬ 1 (2.1)

where {Pi} are geometric coefficients (control points) and the {Bi,n(u)} are the nth-degree
Bernstein polynomials (basis or blending functions) defined as

Bi,n(u)=
n!

i!(n− i)!
ui(1−u)n−i (2.2)



Free vibration and buckling analysis of higher order... 455

In order to accurately represent conic sections, it is necessary to use rational functions instead
of polynomials, so the nth-degree rational Bezier curves are defined as follows

C(u)=

n∑
i=0

Bi,n(u)wiPi

n∑
i=0

Bi,n(u)wi

0¬ u ¬ 1 (2.3)

where by {wi}wemarked the weights that are scalar quantities.
The B-spline curve is a generalization of the Bezier curve defined as

C(u)=
n∑

i=0

Ni,p(u)Pi a ¬ u ¬ b (2.4)

where {Pi} are control points and {Ni,p(u)} are the pth-degree B-spline basis functions (Fig. 1)
that are defined on the nonuniform knot vector

U= [a,. . . ,a
︸ ︷︷ ︸
p+1

,up+1, . . . ,um−p−1,b, . . . ,b︸ ︷︷ ︸
p+1

] (2.5)

The B-spline basis functions of the pth-degree are defined recursively

Ni,0(u)=

{
1 ui ¬ u < ui+1
0 otherwise

Ni,p(u)=
u−ui
ui+1−ui

Ni,p−1(u)+
ui+p+1−u

ui+p+1−ui+1
Ni+1,p−1(u)

(2.6)

Fig. 1. Non-zero linear, quadratic, cubic and quartic B-spline basis functions defined on the open knot

vectorU= [0, . . . ,0︸ ︷︷ ︸
p+1

,0.2,0.5,0.8,1, . . .,1︸ ︷︷ ︸
p+1

]

A rational representation of the B-spline curve is called a NURBS curve. A pth-degree
NURBS curve is defined analogously to (2.3) as

C(u)=

n∑
i=0

Ni,p(u)wiPi

n∑
i=0

Ni,p(u)wi

a ¬ u ¬ b (2.7)

where{Pi} are the control points, {wi} are theweights and {Ni,p(u)} are the pth-degreeB-spline
basis functions that are defined on the nonuniform knot vector given by (2.5).



456 O. Peković et al.

If we define the rational basis functions a

Ri,p(u)=
Ni,p(u)wi
n∑
j=0

Nj,p(u)wj

(2.8)

we can write the NURBS curve as

C(u)=
n∑

i=0

Ri,p(u)Pi (2.9)

It is easy to define multivariate NURBS basis functions by using the tensor product method.
A NURBS surface of degree p in the u direction and degree q in the v direction is a bivariate
vector-valued piecewise rational function of the form

S(u,v) =
n∑

i=0

m∑

j=0

R
p,q
i,j (u,v)Pi,j =

n∑
i=0

m∑
j=0

Ni,p(u)Nj,q(v)wi,jPi,j

n∑
i=0

m∑
j=0

Ni,p(u)Nj,q(v)wi,j

0¬ u,v < 1 (2.10)

where {Pi,j} are the control points, {wi,j} are the weights and {R
p,q
i,j (u,v)} are the bivariate

nonrational B-spline basis function defined on the nonuniform knot vectors

U= [a,. . . ,a
︸ ︷︷ ︸
p+1

,up+1, . . . ,ur−p−1,b, . . . ,b︸ ︷︷ ︸
p+1

] V= [c, . . . ,c
︸ ︷︷ ︸
q+1

,uq+1, . . . ,us−q+1,d, . . . ,d︸ ︷︷ ︸
q+1

]

(2.11)

where r = n+p+1 and s = m+q+1.

Fig. 2. Mesh of control points, the corresponding cubic NURBS curve (left) and NURBS surface (right)

3. Equations of motion

In the TSDT of Reddy (1984, 2004), the displacement field is defined as

u(x,y,z) = u0(x,y)+zψx−
4

3h2
z3
(
ψx+

∂w0
∂x

)

v(x,y,z) = v0(x,y)+zψy −
4

3h2
z3
(
ψy+

∂w0
∂y

)

w(x,y,z) = w0(x,y)

(3.1)

where u0, v0, w0 represent linear displacements of themidplane, ψx, ψy are rotations of normals
to the midplane about the y and x-axes, respectively, and h denotes the total thickness of the
laminate.



Free vibration and buckling analysis of higher order... 457

The in-plane strains {εxx εyy γxy}
T are given as

εp =






εxx
εyy
γxy





=






ε0xx
ε0yy
γ0xy





+z






ε1xx
ε1yy
γ1xy





+z3






ε3xx
ε3yy
γ3xy





=






∂u0
∂x
∂v0
∂y

∂u0
∂y
+
∂v0
∂x






+z






∂ψx
∂x
∂ψy
∂y

∂ψx
∂y
+
∂ψy
∂x






− c1z
3






∂ψx
∂x
+
∂2w0
∂x2

∂ψy
∂y
+
∂2w0
∂y2

∂ψx
∂y
+
∂ψy
∂x
+2

∂2w0
∂x∂y






(3.2)

and the cross plane components {γyz γxz}
T as

εs =

{
γyz
γxz

}

=

{
γ0yz
γ0xz

}

+z2
{
γ0yz
γ0xz

}

=






ψy +
∂w0
∂y

ψx+
∂w0
∂x





−c2z

2






ψy +
∂w0
∂y

ψx+
∂w0
∂x





(3.3)

with c1 =4/(3h
2) and c2 =4/(h

2).
Constitutive relations between stresses and strains in the k-th lamina in the case of plane

stress state, in the local coordinate system of the principle material coordinates (x1,x2,x3),
where x1 is fibre direction, x2 in-plane normal to fibre and x3 normal to the lamina plane, are
given by





σ
(k)
1

σ
(k)
2

τ
(k)
12

τ
(k)
23

τ
(k)
13






=






Q11 Q12 0 0 0
Q12 Q22 0 0 0
0 0 Q66 0 0
0 0 0 Q44 Q45
0 0 0 Q45 Q55






(k)





ε
(k)
1

ε
(k)
2

γ
(k)
12

γ
(k)
23

γ
(k)
13






(3.4)

The quantities Q
(k)
ij are called the plane reduced stiffness components and are given in terms of

material properties of each layer as

Q
(k)
11 =

E
(k)
1

1−ν
(k)
12 ν

(k)
21

Q
(k)
12 =

ν
(k)
12 E

(k)
2

1−ν
(k)
12 ν

(k)
21

Q
(k)
22 =

E
(k)
2

1−ν
(k)
12 ν

(k)
21

Q
(k)
66 = G

(k)
12 Q

(k)
44 = G

(k)
23 Q

(k)
55 = G

(k)
13

(3.5)

E
(k)
1 , E

(k)
2 are Young’s moduli, ν

(k)
12 , ν

(k)
21 are Poisson’s coefficients and G

(k)
12 , G

(k)
23 , G

(k)
13 are shear

moduli of the lamina.
Composite laminates are usually made of several orthotropic layers (laminae) of different

orientation. In order to express constitutive relations in the referent laminate (x,y,z) coordinate
system (Fig. 3), the lamina constitutive relations are transformed as





σ
(k)
xx

σ
(k)
yy

τ
(k)
xy

τ
(k)
yz

τ
(k)
xz






=






Q11 Q12 Q16 0 0
Q12 Q22 Q26 0 0
Q16 Q26 Q66 0 0
0 0 0 Q44 Q45
0 0 0 Q45 Q55






(k)





ε
(k)
xx

ε
(k)
yy

γ
(k)
xy

γ
(k)
yz

γ
(k)
xz






(3.6)



458 O. Peković et al.

where Qij are the lamina plane stress reduced stiffness components in the laminate coordinate
system defined as (Reddy, 2004)






Q11
Q22
Q16
Q26
Q66






(k)

=






m4 n4 2m2n2 4m2n2

n4 m4 2m2n2Q12 4m
2n2

m3 −mn3 mn3−m3n 2(mn3−m3n)
mn3 −m3n m3n−mn3 2(m3n−mn3)
m2n2 m2n2 −2m2n2 (m2−n2)2






(k)




Q11
Q22
Q12
Q66






(k)






Q44
Q45
Q55






(k)

=





m2 n2

−mn mn
n2 m2






(k){
Q44
Q55

}(k)

(3.7)

with m and n denote cosine and sine of the angle θ between the global axis x and local axis x1
(Fig. 3).

Fig. 3. Local and global coordinate systems of a laminate

The dynamic form of the principle of virtual work in matrix form is given by

∫

Ω

δεTpDεp dΩ +

∫

Ω

δεTsD
s
εs dΩ =

∫

Ω

δuTmü dΩ (3.8)

wherem is defined as

m=






I0 0 0 J1 0 −c1I3 0
0 I0 0 0 J1 0 −c1I3
0 0 I0 0 0 0 0
J1 0 0 K2 0 −c1I4 0
0 J1 0 0 K2 0 −c1I4
−c1I3 0 0 −c1I4 0 c

2
1I6 0

0 −c1I3 0 0 −c1I4 0 c
2
1I6






(3.9)

with

(I0,I1,I2,I3,I4,I6)=
N∑

k=1

h/2∫

−h/2

ρ(k)(1,z,z2,z3,z4,z6) dz

J1 = I1− c1I3 K2 = I2−2c1I4+ c
2
1I6

Matrices that relate the stress resultants to the strains are given as

D=





A B E
B D F
E F H




 Ds =

[
As Ds

Ds Fs

]

(3.10)



Free vibration and buckling analysis of higher order... 459

with

(Aij,Bij,Dij,Eij,Fij,Hij)=
N∑

k=1

h/2∫

−h/2

Q
(k)
ij (1,z,z

2,z3,z4,z6)dz

(Asij,D
s
ij,F

s
ij)=

N∑

k=1

h/2∫

−h/2

Q
(k)
ij (1,z

2,z4) dz

uT =

{
u0 v0 w0 ψx ψy

∂w0
∂x

∂w0
∂y

}T

(3.11)

For buckling analysis, the weak form of the virtual work principle in the matrix form is stated
as

∫

Ω

δεTpDεpdΩ +

∫

Ω

δεTsD
s
εs dΩ +h

∫

Ω

[
∇Tδu0 ∇

Tδv0 ∇
Tδw0

]




σ̂0 0 0
0 σ̂0 0
0 0 σ̂0










∇u0
∇v0
∇w0




dΩ

+
h3

12

∫

Ω

[
∇Tδψx ∇

Tδψy
][
σ̂0 0
0 σ̂0

][
∇ψx
∇ψy

]

dΩ =0 (3.12)

where ∇T =
[
∂/∂x ∂/∂y

]
is the gradient operator and σ̂0 =

[
σ0x τ

0
xy

τ0xy σ
0
y

]

is the matrix of

in-plane pre-buckling stresses.

4. Isogeometric finite element formulation of TSDT laminated plate

In TSDT, the field variables are inplane displacements, transverse displacements and rotations
at control points. By using isogeometric paradigm, the same NURBS basis functions that are
used to describe plate geometry are used for the interpolation of field variables

u=






u0
v0
w0
ψx
ψy






=
n×m∑

I=1






NI 0 0 0 0
0 NI 0 0 0
0 0 NI 0 0
0 0 0 NI 0
0 0 0 0 NI











u0I
v0I
w0I
ψxI
ψyI






=
n×m∑

I=1

NIqI (4.1)

where n×m is the number of control points (basis functions), NI are rational basis functions
and qI are degrees of freedom associated with the control point I.
The in-plane strains and shear strains are obtained using Eqs. (3.2),(3.3) and (4.1) as

εp =
∑

I

[B0I +zB
1
I − c1z

3B3I]qI εs =
∑

I

[BS0I − c2z
2BS2I ]qI (4.2)

where

B0 =





N,x 0 0 0 0
0 N,y 0 0 0

N,y N,x 0 0 0




 B1 =





0 0 0 N,x 0
0 0 0 0 N,y
0 0 0 N,y N,x






B3 =





0 0 N,xx N,x 0
0 0 N,yy 0 N,y
0 0 2N,xy N,y N,x






(4.3)



460 O. Peković et al.

and

BS0 =BS2 =

[
0 0 N,y 0 N
0 0 N,x N 0

]

(4.4)

whereN,x andN,y denote thefirst andN,xx,N,yy,N,xy secondderivatives ofN with respect
to x and y.
For free vibration analysis, we can write

(K−ω2M)q=0 (4.5)

and for buckling analysis, we get

(K−λcrKg)q=0 (4.6)

whereK is the global stiffness matrix defined as

K=

∫

Ω





B0

B1

B3






T



A B E
B D F
E F H










B0

B1

B3




+

[
Bs0

Bs2

]T[
As Ds

Ds Fs

][
Bs0

Bs2

]

dΩ (4.7)

The global mass matrixM is given by

M=

∫

Ω

NTmmNm dΩ (4.8)

with

Nm =






NI 0 0 0 0 0 0
0 NI 0 0 0 0 0
0 0 NI 0 0 NI,x NI,y
0 0 0 NI 0 0 0
0 0 0 0 NI 0 0






T

(4.9)

The global geometrical stiffness matrix Kg that takes into account the contribution of shear
strains is given by

Kg =

∫

Ω

NTg IgNg dΩ (4.10)

with

Ng =






∇N 0 0 0 0
0 ∇N 0 0 0
0 0 ∇N 0 0
0 0 0 ∇N 0
0 0 0 0 ∇N
0 0 NI,xx NI,x 0
0 0 NI,xy NI,y 0
0 0 NI,xy 0 NI,x
0 0 NI,yy 0 NI,y






Ig =






Ig0 0 0 0 0 0 0
0 Ig0 0 0 0 0 0
0 0 Ig0 0 0 0 0
0 0 0 Ig2 0 Ig4 0
0 0 0 0 Ig2 0 Ig4
0 0 0 Ig4 0 Ig6 0
0 0 0 0 Ig4 0 Ig6






(4.11)

and

∇N =
[
NI,x NI,y

]T
Ig0 =hσ̂0 Ig2 =

h3

12
σ̂0

Ig4 =
h5

80

[
σ0x −τ

0
xy

−τ0xy σ
0
y

]

Ig6 =
h7

448

[
σ0x −τ

0
xy

−τ0xy σ
0
y

]



Free vibration and buckling analysis of higher order... 461

5. Free vibration analysis of laminated composite plates

In this Section, the performance of quadratic, cubic and quartic TSDT isogeometric elements
in the free vibration analysis is studied. Standard benchmark problems with various plate sha-
pes, boundary conditions, material characteristics and thickness are solved using the proposed
method, and the results are compared to other available ones.

5.1. Square composite plates

We consider [0/90], [0/90/0] and [0/90/0/90] cross-ply laminates. Each ply is made of an or-
thotropic material with the following characteristics: E1 =1.73 ·10

5MPa, E2 =3.31 ·10
4MPa,

G12 = 9.38 · 10
3MPa, G13 = 8.27 ·10

3MPa, G23 = 3.24 ·10
3MPa and ν12 = 0.036. Mass den-

sity ρ is equal to one. The plate is simply supported on all edges, and all layers are assumed
to be of the same thickness and density. The plates length-to-width ratio is a/b = 1 and the
width-to-thickness ratio is b/h =10. The normalized frequency is defined as ω =(ωh)

√
ρ/E2.

InTables 1, 2 and3,wepresent four dimensionless natural frequencies that correspond to the
values of Fourier integers m,n =1,2 for [0/90], [0/90/0] and [0/90/0/90] laminates, respectively.
We compare the results of quadratic, cubic and quartic IGATSDT elements with the exact 3D
elasticity solution and the analytical solutions of TSDT, FSDT and CPT theories given by
Nosier et al. (1993). The shear correction factor for FSDT theory is taken to be π2/12. In this
example, the plate is modeled with 8x8 elements.
The results obtainedwith quadratic elements are in very good agreementwith the analytical

solution based on TSDT theory of Reddy. We see that the results obtained with cubic and
quartic elements are quite similar and differ slightly from the quadratic elements solutions.

Table 1.Non-dimensional frequency parameter ω of the [0/90] laminate

IGA IGA IGA
Exact TSDT FSDT CPT TSDT TSDT TSDT

quadratic cubic quartic

m,n =1 0.06027 0.06057 0.6038 0.06513 0.06056 0.06053 0.06053

m,n =1,2
0.14539 0.14681 0.14545 0.17744 0.14703 0.14664 0.14663

m,n =2,1

m,n =2 0.20229 0.20482 0.20271 0.25814 0.20508 0.20449 0.20448

Table 2.Non-dimensional frequency parameter ω of the [0/90/0] laminate

IGA IGA IGA
Exact TSDT FSDT CPT TSDT TSDT TSDT

quadratic cubic quartic

m,n =1 0.06715 0.06839 0.06931 0.07769 0.06837 0.06835 0.06835

m,n =1,2 0.12811 0.13010 0.12886 0.15185 0.13014 0.12993 0.12993

m,n =2,1 0.17217 0.17921 0.18674 0.26599 0.17957 0.17908 0.17907

m,n =2 0.20798 0.21526 0.22055 0.31077 0.21551 0.21494 0.21494

5.2. Circular composite plates

Next, we consider a symmetric four-layer laminated circular plate with [θ◦/− θ◦/− θ◦/θ◦]
stacking sequence and clamped boundaries.
The material of each ply has the following characteristics: E1 = 40E2, G12 = G13 = 0.6E2,

G23 = 0.5E2, ν12 = 0.25, ρ = 1. Fiber orientation angles are θ = 0
◦,15◦,30◦,45◦ and the



462 O. Peković et al.

Table 3.Non-dimensional frequency parameter ω of the [0/90/0/90] laminate

IGA IGA IGA
Exact TSDT FSDT CPT TSDT TSDT TSDT

quadratic cubic quartic

m,n =1 0.06621 0.06789 0.06791 0.07474 0.06787 0.06785 0.06785

m,n =1,2
0.15194 0.16065 0.16066 0.20737 0.16085 0.16048 0.16048

m,n =2,1

m,n =2 0.20841 0.22108 0.22097 0.29824 0.22133 0.22077 0.22076

diameter-to-thickness ratio is 10. In order to represent the circular plate geometry, we used
quadratic basis functions with knot vectors U = [0,0,0,1,1,1] and V = [0,0,0,1,1,1]. We
chose an appropriate control polygon in order to get a desirable distribution of the parametric
curves. The control polygon and resulting mesh are shown in Fig. 4. The results for a 8× 8
element mesh are presented in Table 4 and compared with the results obtained with MISQ20
elements (Nguyen-Van et al., 2008), MLSDQ elements (Liew et al., 2003) and the IGA FSDT
results obtained by Thai et al. (2012). There is a good agreement between the results. The first
six mode shapes of the quartic [45◦/ − 45◦/ − 45◦/45◦] clamped laminated circular plate are
shown in Fig. 5.

Fig. 4. The control polygon and knot plot of a quadratic circular plate with 4 non-zero knot spans

Fig. 5. First six mode shapes of a quartic [45◦/−45◦/−45◦/45◦] clamped laminated circular plate



Free vibration and buckling analysis of higher order... 463

Table 4.Non-dimensional frequency parameter ω of the [θ◦/−θ◦/−θ◦/θ◦] circular laminated
plate

θ◦ Method
Modes

1 2 3 4 5 6

0◦ FSDT –MISQ20 22.123 29.768 41.726 42.805 50.756 56.950
FSDT –MLSDQ 22.211 29.651 41.101 42.635 50.309 54.553
IGAFSDT quadratic 22.0989 29.5409 40.8126 42.5447 50.2975 54.7732
IGAFSDT cubic 22.1110 29.5550 40.8150 42.5650 50.3201 54.7332
IGAFSDT quartic 22.1227 29.5735 40.8410 42.5854 50.3478 54.7609
IGATSDT quadratic 22.8351 31.4480 45.5659 48.1996 49.5516 56.7189
IGATSDT cubic 22.6745 31.2413 45.3267 48.1985 49.4442 56.5244
IGATSDT quartic 22.6127 31.0741 45.0771 48.1983 49.4310 56.4541

15◦ FSDT –MISQ20 22.698 31.568 43.635 44.318 53.468 60.012
FSDT –MLSDQ 22.774 31.455 43.350 43.469 52.872 57.386
IGAFSDT quadratic 22.6500 31.3012 43.3124 43.3833 52.8952 57.8347
IGAFSDT cubic 22.6626 31.3166 43.3335 43.3899 52.9197 57.8064
IGAFSDT quartic 22.6751 31.3359 43.3550 43.4165 52.9486 57.8349
IGATSDT quadratic 23.4537 33.6251 48.4304 49.4626 58.8442 66.4838
IGATSDT cubic 23.2857 33.4227 48.1945 49.3160 58.5463 66.1343
IGATSDT quartic 23.2140 33.2644 47.9875 49.3001 58.4857 65.9376

30◦ FSDT –MISQ20 24.046 36.399 44.189 52.028 57.478 67.099
FSDT –MLSDQ 24.071 36.153 43.968 51.074 56.315 66.220
IGAFSDT quadratic 23.9428 35.9896 43.7948 50.9574 56.6770 66.0745
IGAFSDT cubic 23.9565 36.0085 43.8164 50.9745 56.7038 66.1011
IGAFSDT quartic 23.9703 36.0298 43.8390 51.0024 56.7337 66.1316
IGATSDT quadratic 24.9036 38.7086 48.9210 56.0703 62.7850 75.2087
IGATSDT cubic 24.7076 38.5058 48.7678 55.8127 62.4374 74.7126
IGATSDT quartic 24.6221 38.4000 48.7367 55.7171 62.3863 74.6206

45◦ FSDT –MISQ20 24.766 39.441 43.817 57.907 57.945 66.297
FSDT –MLSDQ 24.752 39.181 43.607 56.759 56.967 65.571
IGAFSDT quadratic 24.6335 38.9379 43.4120 56.8708 56.9251 65.2751
IGAFSDT cubic 24.6478 38.9591 43.4330 56.8937 56.9531 65.3002
IGAFSDT quartic 24.6622 38.9814 43.4559 56.9205 56.9844 65.3320
IGATSDT quadratic 25.6205 41.4886 48.2065 59.9176 65.6484 73.5627
IGATSDT cubic 25.4140 41.3547 48.0552 59.6995 65.2816 73.0792
IGATSDT quartic 25.3282 41.3150 48.0050 59.6592 65.2714 73.0149

6. Buckling analysis of composite plate

6.1. Square plate under uniaxial compression

We consider a symmetric four-layer [0◦/90◦/0◦/90◦] cross-ply plate with simply supported
(SS-1) boundaryconditions on all sides (Fig. 6).Theplatematerial is the sameas in theprevious
example. In Table 5, we present the convergence study of a dimensionless buckling load factor
defined as λ = λcra

2/(E2h
3)with the edge-to-thickness ratio equal to 10, wherea is edge length,

h is total thickness of the laminate, λcr is the critical load factor andE2 is the elasticmodulus. In
Table 6, the results for different edge-to-thickness ratios and an8×8 elementmesh are compared
with the analytical solutions based on CPT, FSDT and TSDT theories given by Reddy (2004)
andwith IGAFSDT solutions byThai et al. (2012). The obtained results agree remarkablywith
the other available ones.



464 O. Peković et al.

Fig. 6. Simply supported laminated plate under uniaxial (left) and biaxial (right) compression

Table 5.Normalized critical buckling load of the simply supported cross-ply [0◦/90◦/90◦/0◦]

Method
Mesh

4×4 8×8 16×16 32×32

IGATSDT quadratic 23.2336 23.1725 23.1596 23.1563

IGATSDT cubic 23.1563 23.1551 23.1551 23.1551

IGATSDT quartic 23.1551 23.1551 23.1551 23.1551

Table 6. Normalized critical buckling load of the simply supported cross-ply [0◦/90◦/90◦/0◦]
plate

Method
a/h

5 10 20 50 100

CPT (Khdeir and Librescu, 1988) 36.160 36.160 36.160 36.160 36.160

FSDT (Khdeir and Librescu, 1988) 11.575 23.453 31.707 35.356 35.955

TSDT (Khdeir and Librescu, 1988) 11.997 23.340 31.660 35.347 35.953

IGA FSDT quadratic (Thai et al., 2012) – 23.6599 31.8288 35.3945 36.0130

IGA FSDT cubic (Thai et al., 2012) – 23.6594 31.8267 35.3813 35.9617

IGA FSDT quartic (Thai et al., 2012) – 23.6594 31.8267 35.3813 35.9616

IGATSDT quadratic 11.8270 23.1558 31.5738 35.3480 35.9807

IGATSDT cubic 11.8135 23.1386 31.5541 35.3245 35.9474

IGATSDT quartic 11.8134 23.1385 31.5540 35.3243 35.9468

6.2. Square plate under biaxial compression

The last numerical example in this paper considers a symmetric [0◦/90◦/0◦] three-layer sim-
ply supportedplatewith the samematerial characteristics as in the previous example, subjected
to the biaxial buckling load (Fig. 6). the dimensionless buckling factor is defined in the sameway
as in the uniaxial compression example. The results presented inTable 7 show the dimensionless
buckling factor for different length-to-thickness ratios. The results obtained by the proposed
method are in good agreement with CPT, FSDT and TSDT solutions by Khdeir and Librescu
(1988).

7. Conclusions

The current investigation presents the isogeometric free vibration and buckling analysis of a
laminated plate based on the TSDT theory of Reddy. TSDT is chosen in order to avoid the
usage of shear correction factors. By usingNURBSbasis functions, theC1 continuity needed for
the implementation of TSDT inFEM is easily achieved. It is relatively easy and straightforward
to change the order ofNURBSbasis functions, so quadratic, cubic, quartic or higher orderTSDT
elements are easily formulated. The presented results are very accurate and close to analytical



Free vibration and buckling analysis of higher order... 465

Table 7. Normalized critical buckling load of the simply supported [0◦/90◦/0◦] cross-ply plate
under biaxial compression

Method
a/h

5 10 20 50 100

CPT (Khdeir and Librescu, 1988) 14.704 14.704 14.704 14.704 14.704

FSDT (Khdeir and Librescu, 1988) 1.427 5.492 10.202 12.192 13.146

TSDT (Khdeir and Librescu, 1988) 1.465 5.526 10.259 12.226 13.285

IGATSDT quadratic 1.4262 5.2755 9.7590 11.9065 12.9697

IGATSDT cubic 1.4198 5.2670 9.7455 11.8873 12.9437

IGATSDT quartic 1.4198 5.2670 9.7453 11.8868 12.9428

TSDT solutions. It can be seen that for the free vibration and buckling analyses, cubic and
quartic elements give very similar results, while quadratic elements provide slightly less accurate
results. In our opinion, one can use cubic TSDT elements in order to obtain satisfactory results
in the least computationally expensive way.
Isogeometric analysis is not confined only to NURBS basis functions. Since NURBS basis

functions have several disadvantages from the point of view of analysis, such as the inability of
local refinement, a considerable effort is invested in the research of T-splines (Bazilevs et al.,
2010), locally refined (LR) B-splines (Johannessen et al., 2014), PHT splines (Nguyen-Thanh
et al., 2011), hierarchical refinement of NURBS (Schillinger et al., 2012) and other technologies
that are capable of local refinement in the context of isogeometric analysis.
In this paper, the proposed method is used on simple geometries, but it is possible to deal

withmore complex geometries byusingT-spline technique or thebending stripmethodproposed
by Kindl et al. (2010).

Acknowledgment

Thiswork has been supported by theMinistry of Science andTechnologicalDevelopment ofRepublic

of Serbia through Technological Development Project No. 35035.

References

1. Bazilevs Y., Calo V.M.,Cottrell J.A., Evans J.A., Hughes T.J.R., Lipton S., Scott
M.A., Sederberg T.W., Isogeometric analysis using T-splines, Computer Methods in Applied
Mechanics and Engineering, 199, 229-263

2. Benson D.J., Bazilevs Y., Hsu M.C., Hughes T.J.R., 2011, Isogeometric shell analysis:
the Reissner-Mindlin shell, Computer Methods in Applied Mechanics and Engineering, 199, 5/8,
276-289

3. CasanovaC., GallegoA., 2013,NURBS-based analysis of higher-order composite shells,Com-
posite Structures, 104, 125-133

4. Cottrell J.A., HughesT.J.R., BazilevsY., 2009, Isogeometric Analysis: Toward Integration
of CAD and FEA, JohnWiley & Sons, Chichester

5. Echter R., Oesterle B., Bischoff M., 2013, A hierarchic family of isogeometric shell finite
elements,Computer Methods in Applied Mechanics and Engineering, 254, 170-180

6. Hughes T.J.R., Cottrell J.A., Bazilevs Y., 2005, Isogeometric analysis: CAD, finite ele-
ments, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics
and Engineering, 194, 39/41, 4135-4195

7. JohannessenK.A.,KvamsdalT., DokkenT., 2014, Isogeometric analysis usingLRB-splines,
Computer Methods in Applied Mechanics and Engineering, 269, 471-514



466 O. Peković et al.

8. Jones R.M., 1999,Mechanics of Composite Materials, 2nd ed., Taylor & Francis, Philadelphia

9. Khdeir A.A, Librescu L., 1988, Analysis of symmetric cross-ply laminated elastic plates using
a higher-order theory: Part II - Buckling and free vibration,Composite Structures, 9, 4, 259-277

10. Kiendl J., BazilevsY.,HsuM.-C.,WchnerR.,BletzingerK.-U., 2010,The bending strip
method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches,
Computer Methods in Applied Mechanics and Engineering, 199, 37/40, 2403-2416

11. Kiendl J., BletzingerK.-U., Linhard J.,WchnerR., 2009, Isogeometric shell analysiswith
Kirchhoff–Love elements, Computer Methods in Applied Mechanics and Engineering, 198, 49/52,
3902-3914

12. Liew K.M., Huang Y.Q., Reddy J.N., 2003, Vibration analysis of symmetrically laminated
plates based on FSDT using the moving least squares differential quadrature method, Computer
Methods in Applied Mechanics and Engineering, 192, 2203-2222

13. Nguyen-ThanhN.,KiendlJ.,Nguyen-XuanH.,WchnerR.,BletzingerK.U.,Bazilevs
Y.,RabczukT., 2011,Rotation free isogeometric thin shell analysisusingPHT-splines,Computer
Methods in Applied Mechanics and Engineering, 200, 47/48, 3410-3424

14. Nguyen-Van H., Mai-Duy N., Tran-Cong T., 2008, Free vibration analysis of laminated
plate/shell structures based on FSDT with a stabilized nodal-integrated quadrilateral element,
Journal of Sound and Vibration, 313, 1/2, 205-223

15. Nosier A., Kapania R.K., Reddy J.N., 1993, Free vibration analysis of laminated plates using
a layerwise theory,AIAA Journal, 31, 12, 2335-2346

16. Piegl L., Tiller W., 1997,The NURBS Book, Springer-Verlag NewYork, NewYork

17. ReddyJ.N., 1984,Asimplehigher-order theory for laminatedcompositeplates,Journal ofApplied
Mechanics, 51, 4, 745-752

18. Reddy J.N., 2004,Mechanics of Laminated Composite Plates and Shells Theory and Anlysis, 2nd
ed., CRCPress, Boca Raton

19. Rogers D., 2001, An Introduction to NURBS With Historical Perspective, Morgan Kaufmann
Publishers, San Francisco

20. Schillinger D., Dedè L., Scott M.A., Evans J.A., Borden M.J., Rank E., Hughes
T.J.R., 2012, An isogeometric design-through-analysis methodology based on adaptive hierarchi-
cal refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Computer
Methods in Applied Mechanics and Engineering, 249/252, 116-150

21. ShojaeeS.,ValizadehN., IzadpanahE.,BuiT.,VuT.-V., 2012,Freevibrationandbuckling
analysis of laminated composite plates using theNURBS-based isogeometricfinite elementmethod,
Composite Structures, 94, 5, 1677-1693

22. Thai C., Ferreira A.J.M., Carrera E., Nguyen-Xuan H., 2013, Isogeometric analysis of
laminated composite and sandwich plates using a layerwise deformation theory,Composite Struc-
tures, 104, 196-214

23. Thai C.H., Nguyen-Xuan H., Nguyen-Xuan N., Le T-H., Nguyen-Thoi T., Rabczuk
T., 2012, Static, free vibration, and buckling analysis of laminates composite Reissner-Mindlin
plates using NURBS-based isogeometric approach, International Journal for Numerical Methods
in Engineering, 91, 6, 571-603

Manuscript received February 21, 2014; accepted for print November 30, 2014