Isogeometric nonlinear bending and buckling analysis
of variable-thickness composite plate structures
T. Le-Manha , Q. Huynh-Vanb , Thu D. Phanc , Huan D. Phand , H. Nguyen-Xuane,∗

arXiv:1604.01367v1 [] 3 Apr 2016

a

Faculty of Civil Engineering, HCMC University of Technology (HUTECH), Vietnam
b
Faculty of Mechanical Engineering, Nguyen Tat Thanh University, Vietnam
c
Faculty of Mechanical Engineering, Ho Chi Minh City Vocational College, Vietnam
d
Faculty of Mechanical Engineering, Ho Chi Minh City University of Technology, Vietnam
e
Duy Tan university, Danang, Vietnam

Abstract
This paper investigates nonlinear bending and buckling behaviours of composite
plates characterized by a thickness variation. Layer interfaces are described as
functions of inplane coordinates. Top and bottom surfaces of the plate are symmetric about the midplane and the plate could be considered as a flat surface in
analysis along with thickness parameters which vary over the plate. The variable
thickness at a certain position in the midplane is modeled by a set of control points
(or thickness-parameters) through NURBS (Non-Uniform Rational B-Spline) basic functions. The knot parameter space which is referred in modelling geometry
and approximating displacement variables is employed for approximating thickness, simultaneously. The use of quadratic NURBS functions results in C1 continuity of modeling variable thickness and analyzing solutions. Thin to moderately
thick laminates in bound of first-order shear deformation theory (FSDT) are taken
into account. Strain-displacement relations in sense of von-Karman theory are
employed for large deformation. Riks method is used for geometrically nonlinear
analysis. The weak form is approximated numerically by the isogeometric analysis (IGA), which has been found to be a robust, stable and realistic numerical tool.
Numerical results confirm the reliability and capacity of the propose method.
Keywords: Composite plates, Variable thickness, Nonlinear analysis,
Isogeometric analysis

∗

Corresponding author. Email: nguyenxuanhung@duytan.edu.vn

Preprint submitted to Elsevier

February 21, 2018

1. Introduction
Composite materials are extensively employed in industry varying from large
scale structures such as aircraft bodies, ship hulls to specific components of transport vehicles. Particularly, variable thickness panels are widely applicable and
hence requiring appropriate approaches for analysis and design. A number of
studies investigated isotropic beams and plates with the thickness change. Analytical analysis of flexural bars of uniform and variable thickness was presented by
Fertis [1]. Ohga et al. [2] presented the transfer-matrix method for buckling analysis of thin-walled structures with variable thickness cross sections. Timoshenko
and Woinowsky-Krieger [3] gave solutions of circular plates of nonuniform thickness in some particular cases. Wu and Liu [4] investigated free vibration analysis
of variable thickness circular plates using the generalized differential quadrature
rule. Buckling analysis of rectangular plates of linearly tapered thickness was
considered by Eisenberger and colleagues [5, 6] using the extended Kantorovich
method. Theories for analysis of variable thickness composite plates were also
presented. Ganesan and Liu [7] carried out the nonlinear buckling analysis and
the prediction of the first-ply failure of tapered laminated plates under uniaxial compression. Besides, a number of works investigated laminated composite
plates with thickness change caused by ply drop [8, 9]. However, most of studies considered the linear change of thickness and there is a lack of approaches
for composite plates with complex nonuniform thickness. Accordingly, in this
work, a numerical method considering the laminated composite plates of variable
thickness is proposed based on isogeometric analysis in which the thickness of the
plates is approximated from a set of control thickness parameters using NURBS
basis functions.
The idea of using Non-Uniform Rational B-Spline basis functions in finite element analysis procedure was presented by Hughes et al. [10, 11] in so-called isogeometric analysis. The method filled the gap between geometric design based on
NURBS (see [12]) and structural analysis. The smoothness and high continuity in
NURBS-based finite element meshes promise more accurate solutions. Over the
last decade, many efforts have been investigated in developing the method such
as global refinement techniques [10, 13], local refinement techniques [14, 15, 16,
17, 18], quadrature rules [19, 20]. Analysis of various plate structures including
vibration, bending, buckling, postbuckling and nonlinear thermomechanical stability based on IGA was widely investigated in [21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44] and so on.
Nonlinear bending and buckling analysis of variable thickness composite plates
2

are conducted in this work. Kinematics is based on first-order shear deformation theory (FSDT) in considering the global behavior of structures. von-Karman
strain-displacement relations including geometric imperfections which were presented in the works of Le-Manh and Lee [27, 28] in investigating postbuckling of
laminated composited plates are employed. Governing equations are constructed
based on the virtual displacement principle of energy approach. In the isogeometric analysis framework, the 2-dimensional parameter space which is referred in
constructing geometry and finite element model is also the one of thickness approximation. Surfaces of the plate could be arbitrarily defined, in assumption of
being symmetric about the midplane, and then the thickness will be interpolated
from a set of control thickness parameters using NURBS basis. This provides
a flexible solution in modeling smooth nonuniform surfaces while plate theories
are applicable. C1 continuous requirements in modeling variable thickness and
analyzing solutions are achieved by the use of quadratic NURBS elements. Riks
algorithm is employed for geometrically nonlinear analysis. Numerical problems
consisting of isotropic and laminated composite plates in various boundary conditions are carried out. Variable thickness examples are defined and examined in
different amplitudes in order to demonstrate the applicability as well as accuracy
of the present method.
The article is outlined beginning with a general introduction of isogeometric analysis. A brief of theoretical formulations focusing on the adjustment of
kinematics and the nonuniform thickness treatment is presented in Section 3. Numerical examples describing nonlinear bending and buckling behavior of variable
thickness composite plates are carried out in next Section. Some concluding remarks are drawn at last.
2. Isogeometric analysis
A brief recall of NURBS-based geometry and the framework of isogeometric
analysis are presented in this section.
2.1. NURBS basis functions
A knot vector which defines the parameter space is a set of non-decreasing
real numbers,
Ξ = {ξ1 , ξ2 , ..., ξn+p+1 }

3

(1)

where ξi is a knot and interval [ξ1 , ξn+p+1 ] is called a patch. A knot vector is open
if the first and the last knot have the multiplicity p + 1, and uniform if the other
knots ξi are in uniform space.
Assuming Ξ is an open uniform knot vector, B-Spline basis functions of order
p ≥ 0 are constructed using the following recursive formula [12],



1 , if ξi ≤ ξ ≤ ξi+1
0 , otherwise
ξ − ξi
ξi+p+1 − ξ
Ni,p (ξ) =
Ni,p−1 (ξ) +
Ni+1,p−1 (ξ)
ξi+p − ξi
ξi+p+1 − ξi+1
Ni,0 (ξ) =

(2a)
(2b)

B-Spline basis function has important properties such as,
• Non-negative: Ni,p (ξ) ≥ 0, ∀ξ,
• Ni,p (ξ) is locally supported in the interval [ξi , ξi+p+1 ],
P
• Constitute a partition of unity: ni=1 Ni,p (ξ) = 1, ∀ξ.
Considering a 2-dimensional parameter space, the basis functions are derived as,

pq
Rij
(ξ, η) = Ni,p (ξ) Mj,q (η)

(3)

where i = 1, 2, ..., n and j = 1, 2, ..., m. Ni,p (ξ) and Mj,q (η) are the basis of
order p and q in Ξ = {ξ1 , ξ2 , ..., ξn+p+1 } and H = {η1 , η2 , ..., ηm+q+1 } directions,
respectively.
B-Spline surface is a tensor product of a control net of n×m control points
{Bij } ∈ Rd and 2-dimensional B-Spline basis functions,

S (ξ, η) =

n X
m
X

pq
Rij
(ξ, η) Bij (x, y, z)

i=1 j=1

4

(4)

B-Spline surface is C ∞ continuous inside a knot span [ξi , ξi+1 ] and C p−1
continuous at a single knot ξi . NURBS surface is obtained by projecting a
nonrational B-Spline surface S w (ξ, η) from a homogeneous coordinate in Rd+1
space onto physical Rd space, which has a final form similar to Eq. (4) with the
rational NURBS basis defined by

Ni,p (ξ) Mj,q (η) wij
pq
Rij
(ξ, η) = P
n P
m
Ni,p (ξ) Mj,q (η) wij

(5)

i=1 j=1

where wij ≥ 0 is the weight of control point Bijw = (wij xij , wij yij , wij zij , wij )T
in homogeneous space.
2.2. Isogeometric analysis model
Isogeometric analysis proposed by Hughes et al. [10] employs NURBS basis for geometrical modelling and finite element approximations, simultaneously.
More accurate solutions compared with standard finite element analysis of composite structures are usually obtained due to the higher-order continuity in the
NURBS meshes, while using less total degrees-of-freedom (DOFs). In the analysis domain, dependent displacement variables and other geometric information
could be approximated by the following formula,
s=

nC
X

Rc s c

(6)

c=1

where nC is the number of control points, sc are control parameters (i.e.
displacements of control points, initial deviation of control points or control
thickness parameters, etc) and Rc are 2-dimensional NURBS basis functions.
In finite element subdomain, Eq. (6) could be written as,

se =

nCe
X

Rce sec

(7)

c=1

where nCe is the number of control points per element, sec are control parameters
of corresponding element and Rce are NURBS basis functions evaluated with
respect to parameter coordinates which define the element in physical space.
5

3. Theoretical formulation
In this work, first-order shear deformation is employed to consider thin and
moderately thick composite plates. The kinematics of FSDT theory including geometric imperfections is carried out and constitutive relations of composite plates
of variable thickness are presented in the isogeometric analysis framework.
3.1. Kinematics
Assume that midplane is the origin of material coordinate system of a perfectly
flat plate. The displacement field in FSDT is defined as follows [45, 46],

ux (x, y, z) = u(x, y) + zφx (x, y)
uy (x, y, z) = v(x, y) + zφy (x, y)
uz (x, y, z) = w(x, y)

(8a)
(8b)
(8c)

where (ux , uy , uz ) and (u, v, w) are the displacements of a point in body and
displacements of a point on midplane along (x, y, z) directions, respectively. φx
and φy are the rotations of transverse normal of the mid-plane about the y and x
axes, respectively.
When the initial configuration of the plate contains geometric imperfections, i.e.
initial deviations along out-of-plane direction from being perfectly flat, this small
deflection can be described as a function of inplane coordinates,
ūz (x, y, z) = w̄(x, y)

(9)

In a recent work, Le-Manh and Lee [27] applied geometric mode shape
imperfection in investigating postbuckling behavior of rectangular composite
plates. The effect of imperfection magnitude was taken into account. The
deformed geometry of imperfect plate in FSDT could be depicted in Fig. 1. The
total deflection is the summation of w(x, y) and w̄(x, y). Nonlinear
strains-displacement relations in sense of von-Karman theory including
imperfections can be arranged into the inplane strains, curvatures and shear
strains as follows,

6

=

κ=















γ=

∂u
∂x
∂v
∂y
∂u
∂y

+
+
+

∂φx
∂x
∂φy
∂y
∂φx
∂y


1 ∂w 2
∂ w̄
+ ∂w
2  ∂x 
∂x ∂x
2
1 ∂w
∂ w̄
+ ∂w
2 ∂y
∂y ∂y
∂v
∂w
+ ∂w
+ ∂∂xw̄ ∂w
∂x
∂x ∂y
∂y





+

∂w ∂ w̄
∂x ∂y

(10a)









y
+ ∂φ
∂x
+
φy + ∂w
∂y
∂w
φx + ∂x +

(10b)



∂ w̄
∂y
∂ w̄
∂x


(10c)

3.2. Constitutive relations for composite plates of variable thickness
Here, composite plates whose geometry is symmetric about the midplane are
taken into account, and hence their initial configuration still can be modeled as a
flat surface. The thickness of the layer k at a certain position on the midplane is
a function of inplane coordinates h(k) (x, y) which could be approximated from a
(k)
set of control thickness parameters hc using 2-dimensional NURBS basis:
h(k) =

nCh
X

Rc h(k)
c

(11)

c=1

Basically, the number of control points nC and the number of control thickness
parameters nCh could be different, and the parameter spaces for each of them
could be also separated. However, it is more convenient to use a unique
parameter space in the analysis framework. Accordingly, one control thickness
parameter will be assigned at each control point to define a lamina thickness. For
a composite laminate made of n laminae, there will be n × nC control thickness
parameters in total. The thickness of the laminate is derived as follows

h=

n
X
k=1

h(k) =

n X
nC
X

Rc h(k)
c =

k=1 c=1

nC
X
c=1

Rc

n
X

!
hc(k)

(12)

k=1

Due to the change of thickness, the extensional stiffness matrix A, the
extensional-bending coupling matrix B, bending matrix D and transverse shear
stiffness matrix As vary over the plate and can be computed by

7

A=

n
X

_

k=1
n
X

1
B=
2
D=

1
3

_

Q̄(k) z k+1 − z k

k=1
n
X

(13)


(14)

 3

_
_3
Q̄(k) z k+1 − z k

(15)

Q̄

(k)





_2

_2

z k+1 − z k

k=1

As = K s

n
X

_

_

z k+1 − z k
Q̄(k)
s



(16)

k=1

where Q̄(k) and Q̄s are the global plane-stress stiffness and shear stiffness of
_
lamina (k), z k (k = 1, 2, ..., n + 1) are the coordinates of interfaces along normal
direction determined at a specific position on midplane via Gauss quadrature
points in the integral loop and Ks is shear correction factor.
The constitutive relations are then derived as follows,

 


A B 0
 N 
  
M = B D 0  κ




Q
0 0 As
γ

(17)

where N, M, Q is the vectors of inplane force resultants, moment resultants and
transverse shear force resultants, respectively.

8

4. Numerical examples
In this section, numerical problems including bending and buckling analysis of isotropic and composite plates of variable thickness are investigated to illustrate the application of the proposed approach. These plates vary from thin
to moderately thick in bound of FSDT. The finite element meshes of full plates
are constructed using quadratic NURBS elements in nonlinear parameterizations
which lead to regular meshes in physical space. Reduced integration technique
is employed to evaluate the transverse shear terms in Gauss quadrature. Convergence tests are not shown in details for the sake of brevity. For illustration, 6×6
quadratic NURBS elements are employed for analysis. Mesh refinement will be
stated as necessary. Riks method is used for nonlinear analysis with the tolerance
of displacement convergence δ = 10−3 . In buckling analysis, geometric imperfection is applied on the plates in primary mode shape with the normalized amplitude
∆ = w̄c /a = 10−5 , (see Ref. [27]).
The arc length methods for nonlinear analysis have been well known for the
capacity of tracking equilibrium paths through limit point at which the tangent
stiffness matrix KT becomes singular. Among the methods, Riks algorithm is
widely used. Reddy [47] gave an explanation of the scheme and provided a detail
of computational algorithm. The residual vector is a function of both
displacement u and load parameter λ: R (u + w̄, λ) = K (u + w̄) − λF̄; where
w̄ is the imperfection vector, K is the direct stiffness matrix and F̄ is the
preference load vector. In this study, the tolerance of displacement convergence
is δ = 10−3 .
The non-dimensional material properties are defined as: For isotropic plates: E =
3×106 and ν = 0.25, and for composite plate: E1 /E2 = 25, G12 /E2 = 0.5, G23 /E2
= 0.2, G13 = G12 and ν = 0.25.
4.1. Isotropic square plates of linearly tapered thickness along x direction
The model of a square plate with the thickness h linearly decreasing along x
direction is depicted in Fig. 2. Inplane dimension of the plate is a×a. The tapered
ratio α, i.e. the slope of the tapered surface, is defined by
α=

h̄ − hmin
hmax − h̄
=
a
a

9

(18)

where h̄ is the uniform thickness, α = 0, of the plate which has the same volume
with those of tapered thickness. The top surface of the plate could be described
as the following function

z(x) =

h̄
− αx
2

(19)

Moreover, the bottom surface is geometrically symmetric with the top surface
about the midplane. In this example, isotropic square plates are investigated. The
non-dimensional geometry and uniform thickness are assumed as a=10 and
h̄=0.2 (a/h̄ = 50). In bending test, the plate is clamped at AD and subjected to
increasing uniform transverse load, deflection of point M is measured. In
buckling test, all 4 edges of the plate is simply-supported (SS1), uniaxial
compression load is applied in x direction and central transverse displacement at
origin O is measured:
Clamped at AD : u = v = w = φx = φy = 0

(20)

SS1 at AD and BC : w = φy = 0
SS1 at AB and CD : w = φx = 0

(21a)
(21b)

The transverse displacement w, uniform transverse load q and compression load
N are normalized as
w
h̄
qa4
q̄ =
E h̄4
N a2
λ= 2
π D̄
w̄ =

(22)
(23)
(24)

where D̄ is the flexural rigidity of the plate of uniform thickness , i.e. tapered
ratio α = 0.

10

The solutions are given in Fig. 3 and Fig. 4. In bending test, linear responses of
the cantilever plates under uniform transverse load are obtained. It is also
reasonable that the larger tapered thickness ratio α results in smaller deflection.
Fig. 4 presents buckling and postbuckling path of imperfect isotropic square
plates under uniaxial compression load in x direction. The results are in good
agreement with ones of uniform thickness α = 0 given by Sundaresan et al. [48]
and Le-Manh and Lee [27]. It can be seen that when increasing the slope of
tapered surfaces, the critical load of the plate significantly decreases. This
indicates that stability analysis should be taken into account in designing tapered
structures.
4.2. Isotropic square plates of linearly tapered thickness along the diagonal direction
In this problem, isotropic square plates of tapered thickness in two dimensions
are investigated. The thickness h is gradually changes along the diagonal DB of
the plate, Fig. 5. In this case, the tapered ratio α and the top surface of the plate
could be derived as follows
h̄ − hmin
hmax − h̄
√
√
=
(25)
α=
a 2
a 2
h̄
α
α
z(x, y) = − √ x + √ y
(26)
2
2
2
It should be noted the plates with different tapered ratios have the same volume.
Non-dimensional geometry and normalized parameters are similar to the
previous example. In bending test, edge AD and CD are under simply-supported
(SS2) in Eq. (27). The plate is subjected to increasing uniform transverse load
and corresponding transverse displacement is computed at the corner point B. In
buckling test, boundary conditions SS1 in Eq. (21) are applied and the plate is
under biaxial compression loads in x and y directions:
SS2 at AD and DC : u = v = w = 0

(27)

The results of bending test are shown in Fig. 6. It is interesting to observe the
different behavior in linear and nonlinear analysis. When increasing the tapered
ratio α the linear bending deflection decreases, however, when the transverse
displacement is considerably large, higher tapered ratio results in larger nonlinear
bending deflection. Buckling and postbuckling path of the isotropic square plates
under biaxial compression loads are plotted in Fig. 7. The critical buckling load
notably decreases when the slope of tapered surfaces increases.
11

4.3. Symmetric cross-ply (0/90)s laminated composite square plate of linearly tapered thickness along x direction
Considering a cantilever symmetric cross-ply (0/90)s square laminated composite plate of tapered thickness along x direction. The non-dimensional inplane
dimension and uniform thickness of the plate are a=10 and h̄=0.2 (a/h̄ = 50).
Cross section of the laminate in x − z plane is depicted in Fig. 8. The thickness
of each ply is equal together at any position over the plate. Boundary and load
conditions for bending and buckling analysis are similar to example 4.1. Uniform
transverse load q and compression load N are normalized as follows,
qa4
E2 h̄4
N a2
λ=
E2 h̄3

q̄ =

(28)
(29)

The obtained transverse displacement versus loading curves are shown in Fig. 9
for bending analysis and Fig. 10 for buckling analysis. Behavior of the cantilever
symmetric cross-ply laminate has the similar pattern with the clamped isotropic
plate. The deflection of the cantilever composite plate decreases when the
tapered ratio α increases, and the buckling load also drops considerably.
4.4. Symmetric angle-ply (45/-45)s laminated composite square plate of linearly
tapered thickness along the diagonal direction
A symmetric angle-ply (45/-45)s square plate of linearly tapered thickness
along diagonal direction (see Fig. 5) is taken in account here. Boundary and load
conditions are applied same in example 4.2. The thickness of each ply is equal
together at any position over the plate. The inplane dimensions, uniform thickness
and normalized parameters are similar to example 4.3.
The results of bending test are given in Fig. 11. It can be seen that in linear
bending analysis there is an unpredictable development of the deflection with
respect to the increase of tapered ratio α. Therefore, nonlinear analysis should be
considered in designing plates of tapered thickness, especially in cases of
complex boundary conditions. Besides, critical load of the symmetric angle-ply
(45/-45)s laminate under biaxial compression loads notably drops when
increasing the slope of the tapered surface, see Fig. 12.

12

4.5. Composite square plates of thickness as a sine wave
In this example, nonlinear buckling of composite plates including isotropic
square plate and symmetric cross-ply (0/90)s laminated composite square plates of
variable thickness is examined. The top surface of the plates is described as a sine
wave function, Eq. (30), and the bottom surface is geometrically symmetric with
the top surface about the midplane. The cross section of the plates are illustrated
in Fig. 13. Non-dimensional inplane dimension and the uniform thickness are
a=10 and h̄=0.5 (a/h̄=20).
z(x) = αh̄ sin (2π

nx π h̄
+ )+
a 2 2

(30)

where α is the parameter that controls the amplitude of the sine wave and n is to
change the number of wavelength. The plates are subjected to uniaxial
compression load in y direction. When increasing the number of wavelength, the
plate surfaces are more complex and it is required mesh refinement to obtain
convergence solution. In this works, the mesh of 6×6 quadratic NURBS
elements is used for α = 0, 12×12 quadratic NURBS elements for α 6= 0 and
n=1,2 and 3.
The effects of the amplitude parameter α and the number of wavelength n to the
critical load of the plates will be investigated. The obtained bifurcation behavior
in nonlinear buckling analysis is given in Fig. 14 and Fig. 15 for the isotropic
plate, and Fig. 16 and Fig. 17 for the cross-ply laminate. It can be remarked that
α and n slightly affect the buckling and postbuckling path of the isotropic plates
while the buckling load of the cross-ply laminate significantly fluctuates
depended on these parameters. In details, when the wavelength is fixed at n = 1
and increasing the amplitude of sine wave, Fig. 16, the critical load of the
cross-ply laminate dramatically drops in comparison with the laminate of
uniform thickness. On the other hand, when the amplitude of sine wave is fixed
at α = 0.1 and changing the shape the plates by increasing the number of
wavelength, buckling load of the cross-ply laminate decreases considerably
(n=1), but it rises again and approaches a certain value close to the critical load
in case of uniform thickness (n=2,3), see Fig. 17.
5. Conclusions
Nonlinear bending and buckling analysis of symmetric composite plates of
variable thickness based on the framework of isogeometric analysis were pre13

sented in this work. Control thickness parameters were introduced to approximate the thickness which varies over the plates. This allows us model smooth
non-uniform surfaces while using plate theories for analysis. Numerical examples
were provided to illustrate the accuracy and applicability of the method as well
as present some interesting features in nonlinear analysis of isotropic and laminated composite plates of variable thickness. According to the present method
and obtained numerical solutions, some conclusions could be drawn:
(i) The approach successfully obtains nonlinear bending and buckling behavior
of variable thickness composite plates. Numerical problems of isotropic and
laminated composite plates with different kinds of thickness changes were
defined, giving reliable solutions in comparison with other reference methods.
(ii) Nonlinear analysis should be carried out in designing plates of variable
thickness, especially composite laminates whose critical load could considerably
decrease in comparison with ones of uniform thickness having the same volume.

14

References
References
[1] Fertis DG. Nonlinear Mechanics. CRC Press; 1993.
[2] Ohga M, Shigematsu T, Kawaguchi K. Buckling analysis of thin-walled
members with variable thickness. J Struct Eng 1995;121(6):919–24.
[3] Timoshenko S, Woinowsky-Krieger S. Theory of plates and shells. McGrawhill; 1959.
[4] Wu TY, Liu GR. Free vibration analysis of circular plates with variable
thickness by the generalized differential quadrature rule. Thin-Walled Structures 2001;38:7967–80.
[5] Eisenberger M, Alexandrov A. Buckling loads of variable thickness thin
isotropic plates. Thin-Walled Structures 2003;41(9):871–89.
[6] Shufrin I, Eisenberger M. Stability of variable thickness shear deformable platesfirst order and high order analyses. Thin-Walled Structures
2005;43(2):189–207.
[7] Ganesan R, Liu DY. Progressive failure and post-buckling response of
tapered composite plates under uni-axial compression. Compos Struct
2008;82(2):159–76.
[8] Dinardo MT, Lagace PA. Buckling and postbuckling of laminated composite
plates with ply dropoffs. AIAA Journal 1989;27:1392–8.
[9] Varughese B, Mukherjee A. A ply drop-off element for analysis of tapered
laminated composites. Compos Struct 1997;39:123–44.
[10] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: Cad, finite
elements, nurbs, exact geometry and mesh refinement. Comput Methods
Appl Mech Eng 2005;194:4135–95.
[11] Cottrell JA, Hughes TJR, Bazilevs Y. Isogeometric analysis: toward integration of CAD and FEA. John Wiley and Sons; 2009.
[12] Piegl L, Tiller W. The NURBS book (Monographys in Visual Communication), 2nd edition. Springer-Verlag; 1997.
15

[13] Cottrell JA, Hughes TJR, Reali A. Studies of refinement and continuity in isogeometric structural analysis. Comput Methods Appl Mech Eng
2007;196:4160–83.
[14] Bazilevs Y, Calo VM, Cottrel JA, Evans JA, Hughes TJR, Lipson S, et al.
Isogeometric analysis using t-splines. Comput Methods Appl Mech Eng
2010;199:229–63.
[15] Nguyen-Thanh N, Nguyen-Xuan H, Bordas S, Rabczuk T.
Isogeometric analysis using polynomial splines over hierarchical t-meshes
for two-dimensional elastic solids. Comput Methods Appl Mech Eng
2011;200:1892–908.
[16] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wucher R, Bletzinger KU,
Bazilevs Y, et al. Rotation free isogeometric thin shell analysis using phtsplines. Comput Methods Appl Mech Eng 2011;200:3410–24.
[17] Kleiss SK, Juttler B, Zulehner W. Enhancing isogeometric analysis by a
finite element-based local refinement strategy. Comput Methods Appl Mech
Eng 2012;213-216:168–82.
[18] Schillinger D, Dede L, Scott MA, Evans JA, Borden MJ, Rank E, et al. An
isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of nurbs, immersed boundary methods, and t-spline cad
surfaces. Comput Methods Appl Mech Eng 2012;249-252:116–50.
[19] Hughes TJR, Reali A, Sangalli G. Efficient quadrature for nurbs-based isogeometric analysis. Comput Methods Appl Mech Eng 2010;199:301–13.
[20] Auricchio F, Calabro F, Hughes TJR, Reali A, Sangalli G. A simple algorithm for obtaining nearly optimal quadrature rules for nurbs-based isogeometric analysis. Comput Methods Appl Mech Eng 2012;249-252:15–27.
[21] Luu AT, N. I. Kim NI, Lee J. Nurbs-based isogeometric vibration analysis
of generally laminated deep curved beams with variable curvature. Compos
Struct 2015;119:150–65.
[22] Kapoor H, Kapania RK. Geometrically nonlinear nurbs isogeometric
finite element analysis of laminated composite plates. Compos Struct
2012;94:3434–47.

16

[23] Thai CH, Nguyen-Xuan H, Nguyen-Thanh N, Le TH, Nguyen-Thoi T,
Rabczuk T. Static, free vibration, and buckling analysis of laminated composite reissnermindlin plates using nurbs-based isogeometric approach. Int
J Numer Methods Eng 2012;91:571–603.
[24] Thai CH, Ferreira AJM, Carrera E, Nguyen-Xuan H. Isogeometric analysis
of laminated composite and sandwich plates using a layerwise deformation
theory. Compos Struct 2013;104:196–214.
[25] Nguyen-Xuan H, Thai CH, Nguyen-Thoi T. Isogeometric finite element
analysis of composite sandwich plates using a higher order shear deformation theory. Compos Part B 2013;55:558–74.
[26] Nguyen VP, Nguyen-Xuan H. High-order b-splines based finite elements for delamination analysis of laminated composites. Compos Struct
2013;102:261–75.
[27] Le-Manh T, Lee JH. Postbuckling of laminated composite plates using
nurbs-based isogeometric analysis. Compos Struct 2014;109:286–93.
[28] Le-Manh T, Lee JH. Stacking sequence optimization for maximum strengths
of laminated composite plates using genetic algorithm and isogeometric
analysis. Compos Struct 2014;116:357–63.
[29] Guo Y, Nagy AP, Grdal Z. A layerwise theory for laminated composites in
the framework of isogeometric analysis. Compos Struct 2014;107:447–57.
[30] Guo Y, Ruess M, Grdal Z. A contact extended isogeometric layerwise approach for the buckling analysis of delaminated composites. Compos Struct
2014;116:55–66.
[31] Tran LV, Ferreira AJM, Nguyen-Xuan H. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Composites
Part B: Engineering 2013;51:368–83.
[32] Tran LV, Thai CH, Nguyen-Xuan H. An isogeometric finite element formulation for thermal buckling analysis of functionally graded plates. Finite
Elements in Analysis and Design 2013;73:65–76.
[33] Nguyen-Xuan H, Tran LV, Thai CH, Kulasegaram S, Bordas SPA. Isogeometric analysis of functionally graded plates using a refined plate theory.
Composites Part B: Engineering 2014;64:222–34.
17

[34] Hosseini S, Remmers JJ, Verhoosel CV, Borst RD. An isogeometric continuum shell element for non-linear analysis. Comput Methods Appl Mech
Eng 2014;271:1–22.
[35] Kapoor H, Kapania RK, Soni SR. Interlaminar stress calculation in composite and sandwich plates in nurbs isogeometric finite element analysis.
Compos Struct 2013;106:537–48.
[36] Nguyen-Thanh N, Valizadeh N, Nguyen MN, Nguyen-Xuan H, Zhuang X,
Areias P, et al. An extended isogeometric thin shell analysis based on kirchhofflove theory. Comput Methods Appl Mech Eng 2015;284:265–91.
[37] Kiendl J, Bletzinger KU, Linhard J, Wuchner R. Isogeometric shell analysis with kirchhoff-love elements. Comput Methods Appl Mech Eng
2009;198:3902–14.
[38] Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. Isogeometric shell analysis: The reissner-mindlin shell. Comput Methods Appl Mech Eng
2010;199:276–89.
[39] Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. A large deformation, rotation-free, isogeometric shell. Comput Methods Appl Mech Eng
2011;200:1367–13789.
[40] Echter R, Oesterle B, Bischoff M. A hierarchic family of isogeometric shell
finite elements. Comput Methods Appl Mech Eng 2013;254:170–80.
[41] Kiendl J, F. Auricchio L. Beirao da Veiga CL, Reali. A. Isogeometric collocation methods for the reissnermindlin plate problem. Comput Methods
Appl Mech Eng 2015;284:489–507.
[42] Deng X, Korobenko A, Yan J, Bazilevs Y. Isogeometric analysis of continuum damage in rotation-free composite shells. Comput Methods Appl Mech
Eng 2015;284:349–72.
[43] Nguyen LB, Thai CH, Nguyen-Xuan H. A generalized unconstrained theory and isogeometric finite element analysis based on bezier extraction for
laminated composite plates. Eng Comput 2016;32:1–19.
[44] Tran LV, Phung-Van P, Lee J, Wahab M, Nguyen-Xuan H. Isogeometric analysis for nonlinear thermomechanical stability of functionally graded
plates. Compos Struct 2016;140:655–67.
18

[45] Reddy JN. Mechanics of laminated composite plates and shells, 2nd edition.
CRC Press; 2004.
[46] Huu-Tai T, Choi D. A simple first-order shear deformation theory for laminated composite plates. Compos Struct 2013;106:754–63.
[47] Reddy JN. An introduction to nonlinear finite element analysis. Oxford
University Press, New York; 2004.
[48] Sundaresan P, Singh G, Rao GV. Buckling and postbuckling analysis of
moderately thick laminated rectangular plates. Comput Struct 1996;61:79–
86.

19

List of figures

Figure 1: Deformed geometry of imperfect plate in FSDT.

20

Figure 2: Model of a square plate of linearly tapered thickness along x direction.

21

2 0
α=
α=
α=
α=

1 6

0
0 .0 0 5
0 .0 1 0
0 .0 1 5

q

1 2

8

4

0
0

2

4

6

w

8

1 0

M

Figure 3: Deflection versus uniform transverse load of a cantilever isotropic square plate of tapered
thickness along x direction with different tapered ratio α.

1 .6
1 .4
1 .2

- c r

λx / λx

α= 0

1 .0
0 .8
0 .6
α=
α=
α=
α=
α=

0 .4
0 .2

0 , Q 8 F E (S u n d a re s a n e t. a l [5 4 ])
0 , q u a d r a tic N U R B S ( L e & L e e [2 8 ])
0 , q u a d r a tic N U R B S ( p r e s e n t)
0 .0 0 5
0 .0 1 0

0 .0
0 .0

0 .2

0 .4

w

0 .6

0 .8

1 .0

O

Figure 4: Buckling and postbuckling path of a simply-supported isotropic square plate of tapered
thickness along x direction with different tapered ratio α under uniaxial compression load in x
direction.

22

Figure 5: Model of a square plate of linearly tapered thickness along diagonal direction.

23

2 0
L in e a r b e
α=
α=
α=
α=
N o n lin e a
α=
α=
α=
α=

1 6

q

1 2

n d in
0
0 .0 0
0 .0 0
0 .0 1
r b e n
0
0 .0 0
0 .0 0
0 .0 1

g
3 5 4
7 0 7
0 6 1
d in g
3 5 4
7 0 7
0 6 1

8

4

0
0

2

4

6

w

8

1 0

B

Figure 6: Deflection versus uniform transverse load of an isotropic square plate of tapered thickness along diagonal direction with different tapered ratio α.

1 .6
1 .4
1 .2

λ / λc r

α= 0

1 .0
0 .8
0 .6
α= 0
α= 0 .0 0 3 5 4
α= 0 .0 0 7 0 7

0 .4
0 .2
0 .0
0 .0

0 .2

0 .4

w

0 .6

0 .8

1 .0

O

Figure 7: Buckling and postbuckling path of an isotropic square plate of tapered thickness along
diagonal direction with different tapered ratio α under biaxial compression loads in x and y direction.

24

Figure 8: Cross section of a symmetric cross-ply (0/90)s laminated composite plate of linearly
tapered thickness along x direction in x − z plane.

25

2 0

1 6

q

1 2

8

(0 /9 0 )
S

α=
α=
α=
α=
4

0
0 .0 0 5
0 .0 1 0
0 .0 1 5

0
0

2

4

6

w

8

1 0

M

Figure 9: Deflection versus uniform transverse load of a symmetric cross-ply (0/90)s square laminate of tapered thickness along x direction with different tapered ratio α.

1 .6
1 .4
1 .2

- c r

λx / λx

α=0

1 .0
0 .8
0 .6

(0 /9 0 )
0 .4

S

α= 0
α= 0 .0 0 5
α= 0 .0 1 0

0 .2
0 .0
0 .0

0 .2

0 .4

w

0 .6

0 .8

1 .0

O

Figure 10: Buckling and postbuckling path of symmetric cross-ply (0/90)s square laminate of
tapered thickness along x direction with different tapered ratio α under uniaxial compression load
in x direction.

26

1 0
(4 5 /-4 5 )
S

L in e a r b e
α=
α=
α=
α=
N o n lin e a
α=
α=
α=
α=
8

q

6

4

n d in
0
0 .0 0
0 .0 0
0 .0 1
r b e n
0
0 .0 0
0 .0 0
0 .0 1

g
3 5 4
7 0 7
0 6 1
d in g
3 5 4
7 0 7
0 6 1

2

0
0

2

4

6

w

8

1 0

B

Figure 11: Deflection versus uniform transverse load of a symmetric angle-ply (45/-45)s square
laminate of tapered thickness along the diagonal direction with different tapered ratio α.

1 .6
1 .4
1 .2

r

λ / λc

α= 0

1 .0
0 .8
0 .6

(4 5 /-4 5 )
0 .4

S

α= 0
α= 0 .0 0 3 5 4
α= 0 .0 0 7 0 7

0 .2
0 .0
0 .0

0 .2

0 .4

w

0 .6

0 .8

1 .0

O

Figure 12: Buckling and postbuckling path of a symmetric angle-ply (45/-45)s square laminate of
tapered thickness along the diagonal direction with different tapered ratio α under biaxial compression load in x and y directions.

27

Figure 13: Cross section of the plates in x − z plane

28

1 .6
1 .4
1 .2

- c r

λy / λy

α=0

1 .0
0 .8
0 .6

n = 1 : α= 0
n = 1 : α= 0 .1
n = 1 : α= 0 .2

0 .4
0 .2
0 .0
0 .0

0 .2

0 .4

0 .6

w

0 .8

1 .0

O

Figure 14: Buckling and postbuckling path of the isotropic square plate of sine-wave thickness
with fixed the wavelength n = 1 and varying amplitude λ under uniaxial compression load in y
direction.

1 .6
1 .4
1 .2

- c r

λy / λy

α=0

1 .0
0 .8
0 .6
α=
α=
α=
α=

0 .4
0 .2

0
0 .1 : n = 1
0 .1 : n = 2
0 .1 : n = 3

0 .0
0 .0

0 .2

0 .4

0 .6

w

0 .8

1 .0

O

Figure 15: Buckling and postbuckling path of the isotropic square plate of sine-wave thickness
with fixed amplitude λ=0.1 and varying the wavelength n under uniaxial compression load in y
direction.

29

1 .6
1 .4
1 .2

- c r

λy / λy

α=0

1 .0
0 .8
0 .6

(0 /9 0 )

0 .4

S

n = 1 : α= 0
n = 1 : α= 0 .1
n = 1 : α= 0 .2

0 .2
0 .0
0 .0

0 .2

0 .4

0 .6

w

0 .8

1 .0

O

Figure 16: Buckling and postbuckling path of symmetric cross-ply (0/90)s laminated composite
square plate of sine-wave thickness with fixed the wavelength n = 1 and varying amplitude λ
under uniaxial compression load in y direction.
1 .6
1 .4
1 .2

- c r

λy / λy

α=0

1 .0
0 .8
0 .6

(0 /9 0 )
S

α=
α=
α=
α=

0 .4
0 .2

0
0 .1 : n = 1
0 .1 : n = 2
0 .1 : n = 3

0 .0
0 .0

0 .2

0 .4

0 .6

w

0 .8

1 .0

O

Figure 17: Buckling and postbuckling path of the isotropic square plate of sine-wave thickness
with fixed amplitude λ = 0.1 and varying the wavelength n under uniaxial compression load in y
direction.

30

