Composite Structures 106 (2013) 537–548

Contents lists available at SciVerse ScienceDirect

Composite Structures
journal homepage: www.elsevier.com/locate/compstruct

Interlaminar stress calculation in composite and sandwich plates

in NURBS Isogeometric finite element analysis
Hitesh Kapoor a,⇑, Rakesh K. Kapania b, Som R. Soni c
a
Aerospace & Ocean Engineering Department, Virginia Tech, Blacksburg, VA 24060, United States
b
Aerospace & Ocean Engineering Department, 213E Randolph Hall, Virginia Tech, Blacksburg, VA 24061, United States
c
AdTech Systems Research Inc., Dayton, Ohio, United States

a r t i c l e i n f o a b s t r a c t

Article history: This research paper describes the development of NURBS Isogeometric finite element analysis and post-
Available online 10 June 2013 processor for interlaminar stress calculation in composite and sandwich plates. First-order, shear-
deformable laminate composite plate theory is utilized in deriving the governing equations using a var-
Keywords: iational formulation. Linear, quadratic, higher order and k-refined NURBS elements are constructed and
NURBS Isogeometric analysis numerical validation is performed for laminated composite and sandwich plates. Lagrange finite element
Laminated composite plate suffers from higher order stress gradient oscillations due to Gibbs phenomenon and require alternative
First-order shear deformation theory
stress recovery procedures for accurate interlaminar stress calculations, especially interlaminar normal
Interlaminar stress calculation
Direct integration
stress. In this paper, direct post-processing is performed which computes interlaminar shear and normal
k-Refined NURBS higher order element stresses from higher order gradients of NURBS basis in a single step procedure. Interlaminar stresses are
found to be in an excellent agreement with 3D elasticity solution and FSDT along with k-refinement pro-
cedure of NURBS basis is found to compute equivalent or better interlaminar normal stress than higher-
order shear deformation theory.
Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction coarse mesh for the analysis. The idea behind Isogeometric analysis
is to model the geometry exactly which also serves the basis for the
Inspite of extensive use of finite element methods, the barriers solution space i.e. invoking isoparametric concept. Kagan et al. [12]
between engineering design and analysis still exist and the way to developed B-Spline finite element analysis and integrated with
bridge gap is to reconstitute the entire process. Idea is to use one geometric design. Same authors [13] included adaptive refinement
model and use it as an analysis model which require a change from like hp and h-refinement techniques in their finite element code.
classical finite element to an analysis procedure based CAD repre- Advance multi-layered composite and sandwich plate/shell
sentation, formally known as Isogeometric analysis. There are sev- structures are being increasingly used in aerospace, shipbuilding,
eral CAD functions which can be used for CAD representation of bridges and other industries. These structures have smaller thick-
analysis module. Of most widely used CAD basis in engineering de- ness as compared to other dimensions and therefore, are often sub-
sign process are NURBS (non-uniform rational B-splines) as pre- jected to large deformation behavior under external loads. The
sented by Piegle and Tiller [1], Farin [2], Cohen et al. [3] and plate theories for laminated composites can be categorized into
Rogers [4]. NURBS basis are useful for analysis purposes because equivalent single-layer theories (ELS) and layerwise theories.
they possess useful mathematical property of refinement through Equivalent single layer theories namely, classical, first-order and
knot insertion and variational diminishing property of convex hull. higher-order shear deformation theories, are derived from their
There are other computational geometry technologies that can be 3D counterpart (i.e. layerwise and 3D elasticity) by making appro-
utilized as the basis for Isogeometric analysis such as sub-division priate assumptions to the state of strain/stress in the thickness
surface by Peters and Reif [5] and Warren and Weimer [6] Gordon direction, reducing 3D continuum problem to a 2D problem.
patches [7], Gregory patch [8], S-patch [9] and A patch [10], etc. Several articles are available in literature on shear deformation
Hughes et al. [11] introduced the idea of Isogeometric analysis theories. Reissner–Mindlin theory [14–16] (first-order shear defor-
using NURBS (non-uniform rational B-spline). They used NURBS mation theory, FSDT) assumes constant state of through the thick-
to exactly represent the CAD geometry and then, constructed a ness strain. This theory was further extended for anisotropic plates
by Whitney and Pagano [17]. Urthaler and Reddy [18] developed a
⇑ Corresponding author. mixed finite element for bending analysis of laminated composite
E-mail address: hkapoor@vt.edu (H. Kapoor). plates using FSDT. They treated bending moment as a field variable

0263-8223/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.compstruct.2013.05.028
538 H. Kapoor et al. / Composite Structures 106 (2013) 537–548

along with displacement and rotation. Similarly, Cen et al. [19], of transverse stresses. They followed the equilibrium based stress
Kim et al. [20] and Mai-Duy et al. [21] have previously developed recovery method, using the one dimensional, least square finite
elements based on Mindlin–Reissner theory (FSDT). element in the thickness direction.
In higher-order shear deformation theories (HSDT), Putcha and Park and Kim [36] presented a predictor–corrector post-pro-
Reddy [22] developed a mixed finite element approach with 9- cessing procedure for the accurate recovery of stresses and dis-
node Lagrangian quadrilateral element, Kant and Kommineni [23] placement in the multi-layered composite panels. The predictor
formulated refined HSDT with 9-node quadrilateral element for only predicts the transverse shear stress while the corrector meth-
linear and nonlinear finite element analysis of laminated compos- od enhances the accuracy of the displacement, in-plane and trans-
ite and sandwich plates, Polit and Touratier [24] studied large verse normal stress in the thickness direction using the results of
deflection behavior using HSDT, Phan and Reddy [25] developed the predictor and finite element analysis. Park et al. [37], later,
special third order theory (STTR) and Ren and Hinton [26] devel- used the nonlinear predictor–corrector method to obtain the stres-
oped a third order plate theory for the analysis of laminated com- ses and displacement in composite panels in geometrically nonlin-
posite plates. Reddy and Robbins [27]’s provide extensive literature ear formulation. Noor et al. [38] developed a computational
review on ELS and layerwise theories for laminated composite procedure to get the transverse shear stresses in multi-layered
plates. Carrera [28] provide detailed review of zig–zag theories composite panels. They first used the super convergent recovery
for multi layered plates and shells. Comparing higher-order theo- technique to evaluate the in-plane stresses and then, used the
ries with FSDT, higher-order plate theories enhance the accuracy piecewise integration in the thickness direction to obtain the inter-
of the solution slightly but are computationally more expensive. laminar stresses.
On the other hand, FSDT has often been used due to its simplicity Matsunaga [39] analyzed the displacement and stresses in com-
and provide best compromise between economy and accuracy in posite laminated beams using global higher-order beam theory.
predicting the response of thin to moderately thick laminates Author expanded the displacement field variables with power ser-
[29,30]. Kapoor and Kapania [31] developed NURBS Isogeometric ies of z-coordinate. Makeev and Armanios [40] presented an itera-
finite element approach to study geometrically nonlinear behavior tive method to approximate analytical solution of elasticity
of composite plates. In this analysis, they were able to remove problems in composite laminates. Rolfes and Rohwer [41] devel-
locking in thin plates and also were able to formulate a stable oped a method for calculating the improved transverse shear stres-
hourglass behavior using modified NURBS elements in parent ses in laminated composites using first-order shear-deformation
domain. theory. Kant and Manjunatha [42] developed numerical algorithms
It is a well known fact that composites are prone to damage un- for accurate calculations of transverse stresses using higher-order
der low transverse loads due to comparatively low transverse shear-deformation theories. They used direct integration, finite dif-
modulus and have higher probability of failure. The failure or dam- ference and exact surface fitting approach. Recently, Kant et al. [43]
age is generally, initiated due to a variety of failure modes like proposed a semi analytical model for the accurate estimation of
delamination, matrix cracking, fiber failure, etc.; delamination stresses and displacement in composite and sandwich structures.
being the primary mode of failure. Delamination is initiated when The two-point boundary value problem governed by a set of linear
the interlaminar stresses attain the maximum interfacial strength. first-order differential equations through the thickness is solved
Therefore, predicting through-the-thickness stresses accurately is using fourth-order Runge–Kutta method. Nosier and Bahrami
essential. This requires calculating higher-order derivatives of in- [44] developed the analytical solution to study the edge effects
plane stress accurately which in turn requires higher-order dis- in anti-symmetric angle-ply laminates using first order shear
placement derivatives. deformation and Reddy’s layerwise theory.
First order shear deformation theory with the use of constitu- Stress recovery is highly dependent on the structure of a partic-
tive relation produces highly inaccurate interlaminar stresses ular element and on the formulation used in deriving elements.
exhibiting non-physical discontinuity at the ply interface of com- Stress recovery methods include interpolation-extrapolation from
posite laminate. The accurate of interlaminar stresses using 3D super-convergent points [46], L2 projection [47], stress smoothing
equilibrium equations in the context of finite element analysis de- [48] and integral stress techniques [49]. Zienkeiwicz and Zhu
pends upon the computation of higher-order in-plane strain gradi- [50] developed an efficient post-processing technique in terms of
ents, which is generally poor as Lagrange polynomial gradient super-convergent patch recovery (SPR) procedure. A modified ver-
oscillates. Reddy [32] calculated transverse shear stresses using sion of this technique was developed to obtain in-plane stresses at
derivative of in-plane stresses that were obtained by differentiat- nodes and interlaminar stresses using equilibrium equations [52].
ing interpolation functions in finite element approximation. In or- Stress recovery procedures suffer from extraneous stress oscilla-
der to accurately predict delamination, interlaminar normal stress tions [53].
is as important as the shear stress. The computation of interlami- This research present development of NURBS Isogeometric fi-
nar normal stress requires an additional derivative of the basis nite element and interlaminar stress calculations in laminated
function. These higher-order derivatives are not obtained directly composite plates under sinusoidal load. Geometry is modeled ex-
in the finite element code. Tessler [33] developed smoothing vari- actly in Isogeometric framework and isoparametric and sub-para-
ational formulation which combined discrete least square and pen- metric finite element representation is invoked for solution space.
alty constraints functional in a single variational form and First-order, shear-deformable laminated composite plate theory is
recovered the stress gradients more accurately. Byun and Kapania utilized in deriving the governing equations in variational formula-
[34] developed a post-processing technique to overcome this tion. Linear, quadratic, p and k-refined NURBS elements in parent
drawback. They interpolated the finite element displacement data domain are constructed. Interlaminar stresses are computed by di-
using polynomial functions like Chebyshev and orthogonal polyno- rect integration of 3D equilibrium equations, utilizing non-oscilla-
mials in a global domain. Use of Chebyshev polynomials require tory nature of higher order NURBS basis and their gradients. In-
nodal displacement data to be available at some specific points plane stresses are computed using constitutive equations and
i.e. Chebyshev points and cannot interpolate boundary edge nodal strain gradients for interlaminar stress computations are com-
data while orthogonal polynomials use the arbitrarily distributed puted from higher-order gradients of NURBS basis. Most accurate
data points but are more difficult to obtain. Lee and Lee [35] intro- interlaminar stresses, equivalent to 3D elasticity solution, are ob-
duced the non-iterative post processing procedure for the recovery tained by k-refined NURBS elements with coarsest meshes.
 H. Kapoor et al. / Composite Structures 106 (2013) 537–548 539

2. First-order shear deformation plate theory for laminated where N = {Nxx, Nyy, Nxy} are the force resultants, M = {Mxx, Mxy, Myy}
composite plates are the moment resultants and Q = {Qx, Qy} are the shear force resul-
tants. These resultants acting per unit length can be written as:
2.1. Governing equations 8 9 8 9
< Nxx >
> = Z h=2 > < rx >=
The origin of material coordinate system is considered to be the
Nyy ¼ ry dz
>
: >
; h=2 >: >
mid-plane of the laminate and Kirchhoff assumption that the Nxy rxy ;
8 9 8 9
transverse normal remain perpendicular to the mid-surface after < M xx >
> = Z h=2 > < rx > = ð9Þ
deformation is relaxed. However, transverse strain, z, and out-of M yy ¼ ry zdz
plane normal stress, rz are considered negligible. Layers are >
: >
; h=2 >
: >
M xy rxy ;
assumed to be perfectly bonded. The displacement field is   Z h=2  
defined as: Qx sx
¼ dz
Qy h=2 s y
uðx; y; zÞ ¼ u0 ðx; yÞ þ z/x ðx; yÞ
The laminate constitutive equations relate the force and mo-
v ðx; y; zÞ ¼ v 0 ðx; yÞ þ z/y ðx; yÞ ð1Þ
ment resultant to strains in laminate coordinate system through
wðx; y; zÞ ¼ w0 ðx; yÞ ABD matrix.
8 9 2 38 m 9 2 38 b 9
where u0, v0, w0, are the displacement along the x, y and z-axis and < Nxx >
> = A11 A12 A16 >< x >= B11 B12 B16 > < x >=
6 7 6 7
/x, /y are the rotation of transverse normal of the mid-plane about Nyy ¼ 4 A12 A22 A26 5 my þ 4 B12 B22 B26 5 by
>
: >
; >
: m >; >
: b >
the y and x-axis respectively. Nxy A16 A26 A66 xy B16 B26 B66 xy ;
Total in-plane strain  is decomposed into membrane and bend-
ing strains. ð10Þ

 ¼ m þ zb ð2Þ 8 9 2 38 m 9 2 38 b 9
< M xx >
> = B11 B12 B16 <> x >= D11 D12 D16 < > x >=
6 7 m 6 7 b
m, the membrane strain can be written as follows: M yy ¼ 4 B12 B22 B26 5 y þ 4 D12 D22 D26 5 y
8 9 >
: >
; >
: m >; >
: b >
8 m9
> @u0
> M xy B16 B26 B66 D16 D26 D66 xy ;
< x >
> > > xy
= < @x =
@v 0
m
¼  m
y ¼ @y ð3Þ ð11Þ
>
:c >
; >
> >
m
xy
: @u0 þ @ v 0 >
;     
@y @x Qy A44 A45 cyz
¼K ð12Þ
b, the bending strain can be written as follows: Qx A45 A55 cxz
8 9
8 9 > @/y
>
>
< bx >
=
>
>
< @y >
>
= XN 
1   1  
b ¼ > by ¼
@/y
ð4Þ ðAij ; Bij ; Dij Þ ¼ Q kij ðzkþ1  zk Þ; Q kij z2kþ1  z2k ; Q kij z3kþ1  z3k
: >
; >
> >
>
@y
k¼1
2 3
cbxy >
: @/x þ @/y >
;
@y @x N

X 
 ðA44 ; A45 ; A55 Þ ¼ Q k44 ; Q k45 ; Q k55 ðzkþ1  zk Þ ð13Þ
The transverse shear strain can be written as follows: k¼1
  ( @w )
cyz @x
0
þ /x Here, Aij, Bij and Dij represent extensional and shear, extensional-
¼ ð5Þ
cxz @w0
@y
þ /y bending coupling and bending stiffness matrices and Q kij s are
plane-stress reduced stiffness terms. K is the shear correction factor.
The stress strain relation w.r.t material co-ordinate system of a Weak form of the governing equations is obtained by pre-mul-
lamina can be written as: tiplying the equation of motion with du0, dv0, dw0, d /x and d/y
8 9ðkÞ 2 3ðkÞ 8 9
respectively and integrating by parts over the element domain.
< rxx >
> = Q 11 Q 12 Q 16 < xx >
> =
6 7 yy Substituting force, moment and shear force resultant:
ryy ¼ 4 Q 12 Q 22 Q 26 5 ð6Þ
Z  I
>
: >
; >
:c > ;
rxy Q 16 Q 26 Q 66 @du0 @du0
xy 0¼ Nxx þ Nyy dxdy  Px du0 ds
Xe @x @y Ce
 ðkÞ " #ðkÞ   Z  I
syz Q 44 Q 45 cyz @dv 0 @dv 0
¼ ð7Þ 0¼ Nxy þ N yy dxdy  Py du0 ds
sxz Q 45 Q 55 cxz Xe @x @y Ce
Z  
@dw0 @dW 0 @dw0 b @dw0 b @dw0
The equilibrium equations are developed from stress resultants 0¼ Qx þ Q y  dw0 q þ N xx þ N xy
Xe @x @y @x @x @y
by considering balance of force and moment on an infinitesimal  
area of the laminate. The strong form of the governing equation @dw0 b @dw0 b @dw0
þ N xy þ N yy dxdy
of motion can, then, be written as: @y @x @y
I 
@Nxx @N xy b xx @dw0 þ N b xy @dw0 nx
þ ¼0  Qx þ N
@x @y Ce @x @y
@Nxy @N yy  
þ ¼0 b xy @dw 0 b yy @dw 0
@x @y þ Qy þ N þN ny dw0 ds
@x @y
@Q x @Q y Z  I
þ þq¼0 ð8Þ @d/x @d/x
@x @y 0¼ M xx þ Mxy þ d/x Q x dxdy  T x d/x ds
Xe @x @y Ce
@M xx @M xy Z  I
þ  Qx ¼ 0 @d/y @d/y
@x @y 0¼ M xy þ Myy þ d/y Q y dxdy  T y d/y ds
@M xy @M yy Xe @x @y Ce
þ  Qy ¼ 0 ð14Þ
@x @y
540 H. Kapoor et al. / Composite Structures 106 (2013) 537–548

And, the secondary variables are given as: is the knot index where i = 1, 22, . . . , n + p + 1, p is the order of
the polynomial and n is the number of basis functions. The order
P x ¼ Nxx nx þ Nxy ny ; Py ¼ N xy nx þ Nyy ny
of the polynomial, p = 0, 1, 2, 3 . . . , refer to the constant, linear,
T x ¼ M xx nx þ Mxy ny ; T y ¼ Mxy nx þ M yy ny quadratic, cubic piecewise polynomials, respectively. If more than
 
Qn ¼ Qx þ N b xx @w0 þ N
b xy @w0 nx þ Q y þ N b xy @w0 þ N
b yy @w0 ny one knot is located at the same parametric co-ordinate, these are
@x @y @x @y termed as repeated knots. In case of the open knot vector, first
ð15Þ and last knots are repeated p + 1 times. Basis function formed using
the open knot vector are interpolatory at the beginning and end of
where nx and ny are the directional cosines of the outward normal the parametric space interval, i.e. [n1, nn+p+1]. This distinguishes the
vector. knots from the nodes in the finite element analysis as all the nodes
are interpolatory. As a starting point, B-spline basis functions are
3. Linear NURBS Isogeometric finite element formulation defined recursively using Cox–DeBoor algorithm, starting with
piecewise constants (p = 0).
3.1. Displacement field approximation 
1 if ni 6 n < niþ1 ;
Npi ðnÞ ¼ ð21Þ
0 otherwise
The dependent displacement field variables are approximated
by NURBS basis as follows: And, subsequently, basis functions for orders p = 1, 2, 3, . . . , are
X
nCP X
nCP defined as follows:
u0 ðx; yÞ ¼ uj Rj ðxðn; gÞ; yðn; gÞÞ; v 0 ðx; yÞ ¼ v j Rj ðxðn; gÞ; yðn; gÞÞ n  ni p1 niþpþ1  n p1
j¼1 j¼1 Npi ðnÞ ¼ N ðnÞ þ N ðnÞ ð22Þ
niþp  ni i niþpþ1  niþ1 iþ1
X
nCP
w0 ðx; yÞ ¼ wj Rj ðxðn; gÞ; yðn; gÞÞ
j¼1
3.2.2. Rational B-spline (NURBS) basis
X
nCP X
nCP
/x ðx; yÞ ¼ /xj Rj ðxðn; gÞ; yðn; gÞÞ; /y ðx; yÞ ¼ /yj Rj ðxðn; gÞ; yðn; gÞÞ 2D NURBS basis are formed by tensor product of B-spline basis
j¼1 j¼1 in n and g direction and using projective weights associated with
ð16Þ the control points. Rational basis are very similar to B-spline basis
and derive the continuity and support for the function from the
where nCP are the number of control points per element and Rj are knot vector. Also, NURBS basis form a partition of unity and are
2D NURBS basis functions. pointwise non-negative. Rational basis are formed as follows:
Substituting the displacement field approximation into the
weak form, we obtain element stiffness matrix and load vector. Npi ðnÞMqj ðgÞW i;j
Rp;q
i;j ðn; gÞ ¼ PCPw PCPu p
 Rpk fk ¼ 1 : nCPg ð23Þ
For instance, the size of a sub-matrix, K11, in element stiffness i¼1 N
i¼1 i ðnÞM q
j ð gÞW i;j
matrix, is (nCP)  (nCP) and the size of element load vector, F, is
(nCP  dof, 1). nCP stands for control points per element, dof are NURBS surface is defined as,
the no. of degrees of freedom/control point and (a, b = 1, 2, . . . , 5). X
CPw X
CPu
2 11 12 13 14 15
3 Sðx; yÞ ¼ Rp;q
i;j ðn; gÞBi;j ð24Þ
½K  ½K  ½K  ½K  ½K  j¼1 i¼1
6 21 7
6 ½K 
6 ½K 22  ½K 23  ½K 24  ½K 25  7
7
6 7 where Npi ðnÞ and Mqj ðgÞ are the B-Spline basis in n and g direction
½K ab  ¼ 6 ½K 31  ½K 32  ½K 33  ½K 34  ½K 35  7 ð17Þ and Rp;q
6 7 i;j denotes a 2D NURBS basis function. CPw and CPu are the
6 ½K 41  ½K 43  ½K 43  ½K 44  ½K 45  7 number of control points in n and g directions respectively and
4 5
½K 51  ½K 52  ½K 53  ½K 54  ½K 55  nCP = CPu  CPw is the total number of control points per element.
Wi,js are the weights associated with the control points, Bi,j and
Z weights are the vertical coordinates of the corresponding control
F ia¼3 ¼ qRpi dxdy ð18Þ points. Here, weights, Ws considered are equal to one.
Xe

From the element stiffness matrix, the expression for K11, a sub- 3.3. Numerical integration
matrix of element stiffness matrix, K, is as follows:
Z " p p p p
!#
A multi-element patch, a tensor product of knot vectors in n and
11 @Rpi @Rj @Rpi @Rj @Rpi @Rj @Rpi @Rj
K ij ¼ A11 þ A66 þ A16 þ dxdy g directions, i.e. {0 0 0.5 1 1}  {0 0 0.5 1 1}, creates a 2  2 plate
Xe @x @x @y @y @x @y @y @x
mesh in physical domain. NURBS basis multiplied with respective
ð19Þ
control points constructs a geometric model of plate. In order to
The equation of motion can be written in the matrix form as: perform numerical integration, Guass quadrature is employed.
The integration on physical domain is performed by transforma-
½KfDg ¼ fFg ð20Þ
tion from physical to parametric and from parametric to parent do-
main. Mapping between physical and parametric domain is
3.2. NURBS basis performed as follows:
  X
nCP  
x Bxk
This subsection details B-spline and rational basis generation ¼ R1k ðn; gÞ
y k¼1
Byk
and describes the transformation from physical to parametric to
parent domain. The p and k-refined NURBS elements are con- 8 p9 8 @Rp ðn;gÞ 9
structed by modifying the knot vector by knot-insertion. Last part < @Rk = < k =
@x @n
of this subsection describes the construction of various elements. ¼ ½J 1 
: @Rpk ; : @Rpk ðn;gÞ ;
@y @g ð25Þ
" @x @y
#
3.2.1. B-spline basis
@n @n
A knot vector in one dimension is a set of co-ordinates in the ½J ¼ @x @y
parametric space, N = {n1, n2, . . . , nn+p+1}, where ni is the ith knot, i @g @g
 H. Kapoor et al. / Composite Structures 106 (2013) 537–548 541

Bxk and Byk are control point coordinates and Jxnyg is the mapping Similarly, the shape functions in g
^ direction are obtained. The
from physical to parametric domain. The Mapping from parametric product rule yields the following bi-quadratic shape functions:
to parent domain is the standard mapping as is done in finite ele-
ment. Fig. 1 shows schematic of Isogeometric framework. R21 ð^n; g
^ Þ ¼ M 21 ð^nÞN21 ðg
^ Þ R22 ð^n; g
^ Þ ¼ M 21 ð^nÞN22 ðg
^Þ

3.4. NURBS parent elements R23 ð^n; g

^ Þ ¼ M 21 ð^nÞN23 ðg
^Þ

Various lower and higher order NURBS elements are con- R24 ð^n; g
^ Þ ¼ M 22 ð^nÞN21 ðg
^ Þ R25 ð^n; g
^ Þ ¼ M 22 ð^nÞN22 ðg
^Þ
structed by utilizing various refinement techniques. Linear NURBS
R26 ð^n; g
^ Þ ¼ M 22 ð^nÞN23 ðg
^Þ
element, 9LinNURBS1,2,3/F with 9 control point, constructed by a
knot insertion at {0}, {0.25} and {0.5} locations respectively, qua-
dratic NURBS element, 9QuadNURBS, with 9 control points, k-re- R27 ð^n; g
^ Þ ¼ M 23 ð^nÞN21 ðg
^ Þ R28 ð^n; g
^ Þ ¼ M 23 ð^nÞN22 ðg
^Þ
fined quadratic NURBS element, 16QuadNURBSKR and cubic R29 ð^n; g
^ Þ ¼ M 23 ð^nÞN23 ðg
^Þ
NURBS element, 16CubicNURBS, with 16 control point, k-refined
cubic NURBS element, 25CubicNURBSKR, and quartic NURBS ele- Fig. 2 shows NURBS basis in ^
n and g
^ directions along with the parent
ment, 25QuarticNURBS, with 25 control points and k-refined quar- element.
tic NURBS element, 36QuarticNURBSKR, are developed for
interlaminar stress calculations. In element notation, number in 3.4.2. k-refined quadratic NURBS element (16QuadNURBSKR)
front denotes the number of control points/element and KR de- A tensor product of knot vector {1 1 1 0 1 1 1} in ^
n and g
^
notes k-refinement procedure. A few NURBS elements construction direction results in a bi-quadratic NURBS element, 16Quad-
in parent domain are detailed here. NURBSKR1(F/R). The knot vector is divided into two intervals
and shape functions have unique value over each interval. The
3.4.1. Quadratic NURBS element (9QuadNURBS)) shape functions in ^n direction over each interval are derived as
A tensor product of knot vector {1 1 1 1 1 1} in ^ n and g^ follows:
direction results in a bi-quadratic NURBS element in the parent do- ! !
main. The element consists of 9 control points. The shape function ^n4  ^n ^n4  ^n
in ^
n direction can be written as follows: M 21 ð^nÞ ¼
^n4  ^n2 ^n4  ^n3
! !
^n4  ^n ^n4  ^n ! ! ! !
M21 ð^nÞ ¼ ^n  ^n2 ^n4  ^n ^n5  ^n ^n  ^n3
^n4  ^n2 ^n4  ^n3 M 22 ð^nÞ ¼ þ ^n3 6 ^n < ^n4
! ! ! ! ^n4  ^n2 ^n4  ^n3 ^n5  ^n3 ^n4  ^n3
n^  n^2 n^4  n^ ^n5  n^ ^n  ^n3
! !
M22 ð^nÞ ¼ þ ^n3 6 ^n < ^n4
^n  ^n3 ^n  ^n
^n4  ^n2 ^n4  ^n3 ^n5  ^n3 ^n4  ^n3
! ! M 23 ð^nÞ ¼
n^  n^3 ^n  n^ n^5  n^3 ^
^n4  n3
M23 ð^nÞ ¼
^n5  ^n3 ^n4  ^n3 M 24 ð^nÞ ¼ 0
ð26Þ ð27Þ

Fig. 1. Mapping between physical and parent domain: a framework for Isogeometric finite element analysis.
542 H. Kapoor et al. / Composite Structures 106 (2013) 537–548

Fig. 2. Quadratic NURBS element with basis functions in each direction.

M 21 ð^nÞ ¼ 0 4. Calculation of interlaminar stresses

! !
^n5  ^n ^n5  ^n The evaluation of transverse stresses, sxz, syz and rzz from the
M 22 ð^nÞ ¼
^n5  ^n3 ^n5  ^n4 stress–strain constitutive relations leads to discontinuity at the
! ! ! ! interface of the two adjacent layers of the laminate. Three-dimen-
^n  ^n3 ^n  ^n ^n  ^n3 ^n5  ^n sional analysis is very complex due to variation in constitutive laws
M 23 ð^nÞ ¼ þ ^n4 6 ^n < ^n5
n^5  n^3 ^n4  n3^ ^n5  ^n3 ^n5  n^3 in the thickness direction and continuity requirement of transverse
! ! stresses and displacement across the interface. Therefore, 3D elas-
^n  ^n4 ^n4  ^n ticity equilibrium equations are used to calculate the transverse
M 24 ð^nÞ ¼ stresses in the kth lamina and the equations are given as:
n^6  ^n4 n^5  n^4
ð28Þ @ rkxx @ rkxy @ skxz
þ þ ¼0
Similarly, the shape functions in g^ direction are obtained. The @x @y @z
product rule yields 16 bi-quadratic shape functions, R2k , requiring @ rxyx @ ryy @ syz
k k k
þ þ ¼0 ð29Þ
4  4 Gauss points. In k-refinement process, the order of the exist- @x @y @z
ing knot vector from base model is elevated first and then, a knot is
@ skxz @ skyz @ rkzz
inserted instead of elevating the order after knot insertion. This re- þ þ ¼0
@x @y @z
sults in q  1 continuous derivatives instead of p  1 continuous
derivatives, where p and q are the order of the knot vector before The stresses are calculated using the direct integration method.
and after order elevation, respectively. It can be seen that the equations for transverse shear stress, sxz and
For interlaminar stress calculation, two other modifications of syz are a first-order equation and solving these equations result in
k-refined quadratic NURBS elements are formulated by a single one constant of integration. However, this problem is complax as
knot insertion in the knot vector at different locations. 16Quad- the transverse shear stresses are known at the top and bottom of
NURBSKR2 and 16QuudNUUBSKR3 elements are formed by a tensor the laminate. Therefore, one can only obtain a non-unique solution
product of similar knot vector, {1 1 1 0.25 1 1 1} and as both the traction boundary conditions cannot be enforced
{1 1 1 0.5 1 1 1} in ^
n and g^ directions, respectively. Fig. 3 shows simultaneously. For transverse normal stress, a second-order equa-
k-refined, quadratic NURBS element in the parent domain and tion is derived and require double integration through the thick-
NURBS basis in each direction. ness. The transverse stresses are computed directly as follows,
 H. Kapoor et al. / Composite Structures 106 (2013) 537–548 543

Z z
@ rxx @ rxy 5.1. Simply supported, four layer, cross-ply (0/90/90/0) square plate
sxz ¼  þ dz þ C 1
z @x @y
Z kz  In this analysis, a simply-supported, square, cross-ply (0/90/90/
@ rxy @ ryy
syz ¼  þ dz þ C 2 ð30Þ 0) laminated composite plate with a/h = 10 ratio, under sinusoidal
zk @x @y
Z z Z ! ! loading condition is studied. Material for each lamina is considered
@ 2 rxx @ 2 ryy @ 2 sxy as homogeneous, elastic and orthotropic. Material properties used
rzz ¼ þ þ2 dz dz þ zC 3 þ C 4
zk z @x2 @y2 @x@y are as follows: E1 = 25E2, G12 = G13 = 0.5E2, G13 = 0.2E2, m12 = 0.25,
and Ks = 5/6. 9QuadNURBS with 3  3 mesh and 2  2 Gauss inte-
C1, C2, C3 and C4 are the constants of integration. The constants, C1
gration scheme and family of k-refined quadratic NURBS elements,
and C2 are obtained from the known values of sxz and syz at either
16QuadNURBSKR, with 2  2 mesh and 4  4 Gauss integration
of the two boundaries at z = ± h/2 while C3 and C4 are obtained from
scheme are used for computing in-plane and interlaminar stresses.
the known values of rz at z = ± h/2. Also, noteeworthy is that the
The in-plane stresses are obtained using constitutive equations
calculation of rz requires the use of third derivative of
with smoother displacement gradients. Interlaminar stresses are
displacement.
computed from 3D equilibrium equations by direct integration,
without any recovery procedure and are computed from strain
5. Numerical testing gradients.
The computed stresses are compared with 3D elasticity solu-
Different boundary conditions, plate to thickness ratios and tionPagano and 4Q8Reddy/FE i.e. 4  4 mesh of 8 node quadratic La-
symmetric and anti-symmetric, cross-ply and angle ply laminated grange element. Interlaminar stresses are found to be in close
composite and sandwich plate are studied for numerical valida- proximity to 3D elasticity solution for the k-refined quadratic
tion. Due to bi-axial symmetry, only a quadrant of the plate is mod- NURBS elements. The k-refinement provide smoother in-plane
eled and interlaminar stresses are computed at nearest Gauss stress gradients and thus, resulting in smoother interlaminar stres-
points of specific location. The boundary conditions are as follows: ses computation for the coarsest mesh. Figs. 4 and 5 shows com-
SS1 : v 0 ¼ w0 ¼ /y ¼ 0; at x ¼ a=2; u0 ¼ w0 ¼ /x ¼ 0; at y ¼ b=2 parison of interlaminar shear stresses s xz and syz for family of
NURBS elements with 3D elasticity solution. Table 1 shows the
SS2 : u0 ¼ w0 ¼ /y ¼ 0; at x ¼ a=2; v 0 ¼ w0 ¼ /x ¼ 0; at y ¼ b=2
comparison of in-plane and interlaminar stresses computed using
symmetry B:C:; x ¼ 0; u0 ¼ /x ¼ 0 y ¼ 0; v 0 ¼ /y ¼ 0 NURBS post-processor with 3D elasticity solution.

Fig. 3. k-refined Quadratic NURBS element with basis functions in each direction.
544 H. Kapoor et al. / Composite Structures 106 (2013) 537–548

Ks = 5/6. A family of k-refined quadratic NURBS elements, 16Quad-

NURBSKR with 2  2 mesh and two Gauss integration schemes i.e.
3  3 and 4  4, are used.
Interlaminar stresses are computed from direct integration of
3D equilibrium equations, without any recovery procedure. Com-
puted stresses are compared with 3D elasticity solutionPagano, Byun
and Kapania (7  7 mesh of 9 node quadratic Lagrange element
with stress recovery procedure) and 4Q8Reddy/FE (4  4 mesh of
8 node quadratic Lagrange element).
Stresses are found to be in an excellent agreement with 3D elas-
ticity solution for at least one of the elements in the family of k-re-
fined quadratic NURBS elements when compared with 3D
Fig. 4. Interlaminar shear stress r
 xz (a/2, 0, z) in, a simply supported (SS1), cross- elasticity solution and Byun and Kapania’s results. However, to ob-
ply (0/90/90/0) square plate under sinusoidal loading.
tain the most accurate interlaminar stresses, a single knot in the
knot vector of a particular element of the family of k-refined qua-
dratic NURBS elements, for which the most accurate interlaminar
stress is obtained, is replaced by a new knot. This procedure pro-
duces the most accurate interlaminar stresses which are almost ex-
act to 3D elasticity solution. This new knot produces more
smoother in-plane stress gradients and thus, interlaminar stresses.
This trend is seen for all the length to thickness ratios and for both
integration schemes. Tables 2–4 provide the comparison of in-
plane and interlaminar stresses with 3D elasticity solution, Byun
and Kapania’s and Reddy’s results.

5.3. A simply-supported, two layer, cross-ply (0/90), square laminate

Fig. 5. Interlaminar shear stress r

 yz ð0; b=2; zÞ in a simply supported (SS1), cross- A simply supported, two layered (0/90), square laminate with a/
ply (0/90/90/0) square plate under sinusoidal loading. h ratio = 4, subjected to sinusoidal loading is considered. The mate-
rial properties used here are as follows:

5.2. A simply-supported, three layer, cross-ply (0/90/0) square plate E1 =E2 ¼ 25; G12 =E2 ¼ 0:50; G23 =E2 ¼ 0:20; G13 ¼ G12 ; E2
¼ E3 ; m12 ¼ m23 ¼ m13 ¼ 0:25 ð31Þ
In this analysis, a simply-supported, square, cross-ply (0/90/0)
laminated composite plate with various length to thickness ratios, Computed interlaminar stresses are compared with Kant’s and 3D
a/h = 10, 20, 100, subjected to sinusoidal loading is studied. Mate- elasticity solution. Kant uses various techniques such as finite dif-
rial for each lamina is considered as homogeneous, elastic and ference and exact surface fitting to compute interlaminar stresses
orthotropic. Material properties used are as follows: E1 = from in-plane stresses along with HSDT. Cubic and k-refined qua-
1.72369  1011 N/m2, E2 = 6.894757  109 N/m2, G12 = G13 = dratic NURBS elements are developed to compute interlaminar
3.447378  109 N/m2, G23 = 1.378951  109 N/m2, m12 = 0.25, and shear stresses. In k-refinement procedure, a single knot is inserted

Table 1
 of a simply supported (SS1), cross-ply (0/90/90/0) square plate under sinusoidal loading with 3D elasticity solutionPagano.
Comparison of non-dimensionalized stresses r

a/h Source r ð0;0;h=2Þ r ð0;0;h=4Þ ða=2;b=2;h=2Þ

r xy r ða=2;0;0Þ
xz
ð0;b=2;0Þ
r yz
xx yy

10 ELS/Pagano 0.5590 0.4010 0.0276 0.301 0.196

9QuadReddyFE 0.4954 0.3589 0.0240 0.318 0.181
9QuadNURBS 0.5197 0.4029⁄ 0.0263 0.2980⁄ 0.1812
16QuadNURBSKR1 0.5703 0.3959 0.0307 0.4211 0.2606
16QuadNURBSKR2 0.5365 0.4142 0.0283 0.3194 0.2127
16QuadNURBSKR3 0.5478⁄ 0.4193 0.0274⁄ 0.2976⁄ 0.2057⁄

Table 2
Non-dimensionalized stress r
 in a simply supported (SS1), cross-ply (0/90/0) square plate under sinusoidal loading for a/h = 10.

a/h Source r xx r yy r xy r xz r yz
10 16QuadNURBSKR1 0.4658 0.2466 0.0237 0.5109 0.1564 44
16QuadNURBSKR1 0.5928⁄ 0.3118 0.0319 0.1755 0.0559 33
16QuadNURBSKR2 0.5571 0.2859 0.0295 0.3995 0.1345⁄ 44
16QuadNURBSKR2 0.5294 0.2738 0.0257 0.3204⁄ 0.1107⁄ 33
16QuadNURBSKR3 0.5664 0.2873⁄ 0.0287⁄ 0.3587⁄ 0.1378 44
16QuadNURBSKR3 0.4816 0.2470 0.0227 0.7397 0.2348 33
16QuadNURBSKROP – – – 0.3571(0.4985) 0.1246(0.275) 44
16QuadNURBSKROP – – – 0.3575(0.29) 0.1259(0.30) 33
ElasticityPagano 0.590 0.288 0.0290 0.357 0.1228
Byun/Kapania 0.513 0.254 0.0252 0.380 0.1106
Reddy/FE 0.5098 0.2518 0.0250 0.4060 0.0908
 H. Kapoor et al. / Composite Structures 106 (2013) 537–548 545

Table 3
Non-dimensionalized stress r
 in a simply supported (SS1), cross-ply (0/90/0) square plate under sinusoidal loading for a/h = 20.

a/h Source r xx r yy r xy r xz r yz
1
20 16QuadNURBSKR 0.4883 0.1959 0.0221 0.5246 0.1164 44
16QuadNURBSKR1 0.6179 0.2390 0.0263 0.1764 0.0443 33
16QuadNURBSKR2 0.5786 0.2156 0.0250 0.4191 0.1160⁄ 44
16QuadNURBSKR2 0.5513⁄ 0.2088⁄ 0.0231⁄ 0.3305⁄ 0.0880⁄ 33
16QuadNURBSKR3 0.5006 0.2166 0.0253 0.3839⁄ 0.1333 44
16QuadNURBSKR3 0.5890 0.1864 0.0229 0.7641 0.1768 33
16QuadNURBSKROP – – – 0.3851(0.501) 0.0935(0.31) 44
16QuadNURBSKROP – – – 0.3847(0.304) 0.0935(0.275) 33
ElasticityPagano 0.552 0.210 0.0234 0.385 0.0938
Byun and Kapania 0.532 0.200 0.0223 0.390 0.0899
Reddy/FE 0.5281 0.198 0.0222 0.4176 0.0754

Table 4
Non-dimensionalized stress r
 in a simply supported (SS1), cross-ply (0/90/0) square plate under sinusoidal loading for a/h = 100.

a/h Source r xx r yy r xy r xz r yz
100 16QuadNURBSKR1 0.5371⁄ 0.1903 0.0217⁄ 0.4463⁄ 0.0886⁄ 44
16QuadNURBSKR1 0.6256 0.2118 0.0224 0.1750 0.0415 33
16QuadNURBSKR2 0.55668 0.1821⁄ 0.0232 0.5722 0.1459 44
16QuadNURBSKR2 0.5594 0.1853 0.0225 0.3346⁄ 0.0781⁄ 33
16QuadNURBSKR3 0.5788 0.1846 0.0246 0.6651 0.2055 44
16QuadNURBSKR3 0.5069 0.1647 0.0242 0.7754 0.1496 33
16QuadNURBSKROP – 0.181(0.268) – 0.3957(0.0407) 0.0827(0.02) 44
16QuadNURBSKROP – – – 0.3954(0.309) 0.0826(0.275) 33
ElasticityPagano 0.539 0.181 0.0213 0.395 0.0828
Chansup Byun 0.539 0.181 0.0213 0.394 0.0825
Reddy/FE 0.5346 0.1791 0.0212 0.4215 0.0699

Table 5 r yz is found to be within 1% of 3D elasticity solution. In k-refine-

Non-dimensionalized interlaminar shear stresses r  xz;yz in a simply-supported (SS1), ment procedure, a knot insertion at a location 0.515 in the knot
cross-ply (0/90) square plate under sinusoidal loading for a/h = 4. vector of k-refined quadratic NURBS element computes value of
a/h Source r xz r yz Mesh r yz as 0.3161 when compared to 0.3188 obtained from 3D elastic-
1
ity’s solution, 0.319233 obtained from Kant’s HOST11(16N) (direct
4 16QuadNURBSKR – – 22
16QuadNURBSKR2 – – 22 integration (DI)) element, which is a cubic, 16 node element with
16QuadNURBSKR3 – 0.3541 22 11° of freedom per node and 0.32580 computed from Ren’s higher
16QuadNURBSKROP 0.3157(0.515) 0.3161(0.515) 22 order shear deformation theory model. Fig. 6 shows the compari-
16CubicNURBS 0.3403 0.3688 22 son of through the thickness distribution of r  xz with 3D elasticity
16QuadNURBSKROP 0.3137(0.515) 0.3151(0.515) 33
solution. Similarly, r  yz is computed using cubic and quadratic
16CubicNURBS 0.3316 0.3489 33
ElasticityPagano 0.3127 0.3188 NURBS element. Fig. 7 shows the comparison of interlaminar shear
Ren 0.32450 0.32580 stress, r  yz , with 3D elasticity solution.
HOST7A (9N)(FD) 0.313515 0.294540 However, r  zz requires calculation of third-order displacement
HOST12(16N)(FD) 0.316315 0.319463
derivatives and thus, require higher order polynomial approxima-
tion for derivative calculations. Cubic and quartic NURBS elements
between 1 and 1 and the analysis is run for different knot inser- are constructed for computing r  zz . From Table 6, it is observed that
tions to obtain the most accurate stresses. Tables 5 and 6 compares NURBS elements produce better estimate of through thickness dis-
the interlaminar shear and normal stresses with 3D elasticity, tribution of r  zz as compared to Kant’s solution. k-refinement proce-
Kant’s and Ren’s solutions. dure enhances the smoothness of strain gradients and thus, 3D
From Table 6, for cubic NURBS element, the computed r  yz is elasticity equivalent interlaminar normal stresses are obtained
found to be within 9% and for k-refined NURBS element, computed for various single knot insertion values. Fig. 8 shows the

Table 6
Non-dimensionalized interlaminar normal stress r
 zz in a simply-supported (SS1), cross-ply (0/90) square plate under sinusoidal loading for a/h = 4.

z/h 3D Elasticity HOST11(16 N)(DI) 16CubicNURBS 25CubicNURBSKR 25QuarticNURBS 36QuarticNURBSKR

0.5 1 1 1 1 1 1
0.4 0.973684 0.924730 0.9011 0.9725(0.775,0.775) 0.9051 0.9233
0.3 0.960526 0.817625 0.8023 0.9646(0.7765,0.7765) 0.8102 0.84659
0.2 0.921053 0.698580 0.7034 0.9260(0.7755,0.7755) 0.7152 0.76988
0.1 0.868421 0.58460 0.6046 0.8654(0.774,0.774) 0.6204 0.6932
0.0 0.789474 0.491027 0.5057 0.78832(0.7725,0.7725) 0.5255 0.6165
0.1 0.671050 0.399094 0.3954 0.67387(0.75,0.175) 0.3796 0.4427
0.2 0.460526 0.287604 0.2966 0.5053(0.75,0.175) 0.2847 0.3321
0.3 0.25 0.171356 0.1977 0.2568(0.75,0.045) 0.1898 0.2241
0.4 0.07899 0.072587 0.0989 0.07815(0.768,0.75) 0.0949 0.0767
0.5 0.0 0.0 0.0 0.0 0.0 0.0
546 H. Kapoor et al. / Composite Structures 106 (2013) 537–548

Fig. 6. Interlaminar shear stress, r

 xz ða=2; 0; zÞ, in a simply-supported (SS1), cross- Fig. 9. Through-the-thickness distribution of interlaminar shear stress,
ply (0/90) square plate under sinusoidal loading. r xz ða=2; 0; zÞ, in a simply-supported (SS2), anti-symmetric (45/45) square
laminated plate, a/h = 4.

Fig. 7. Interlaminar shear stress r

 yz ð0; b=2; zÞ in a simply-supported (SS1), cross-
ply (0/90) square plate under sinusoidal loading. Fig. 10. Through-the-thickness distribution of interlaminar shea stress, r
 zz ð0; 0; zÞ,
in a simply-supported (SS2), antisymmetric (45/45) square laminated plate, a/
h = 4.

Fig. 8. Interlaminar normal stress r

 zz ð0; 0; zÞ in a simply-supported (SS1), cross-ply
(0/90) square plate under sinusoidal loading.
Fig. 11. Through-the-thickness distribution of interlaminar normal stress,
r zz ð0; 0; zÞ, in a simply-supported (SS2), antisymmetric (45/45) square laminated
plate, a/h = 10.

Table 7
Non-dimensionalized displacement and in-plane stresses in a simply-supported
(SS2), cross-ply (45/45) square plate under sinusoidal loading.

a/h Source   102

w   102
w r xx r xy
(n = 2) (n = 10) (n = 10) (n = 10)
4 9QuadNURBS 1.1850 0.8645 0.14377 0.14016
9LinNURBS 1.18038 0.8654 0.1386 0.1349
Closed Form 1.1576 0.8450 0.1430 0.1392
10 9QuadNURBS 0.5907 0.27098 0.14377 0.1402
9LinNURBS 0.58395 0.2690 0.1386 0.1350
Closed Form 0.5773 0.2647 0.1430 0.1392
20 9QuadNURBS 0.5053 0.1860 0.14379 0.1403
9LinNURBS 0.4987 0.1838 0.13865 0.1350
Closed Form 0.4944 0.1818 0.1430 0.1392
100 9QuadNURBS 0.4757 0.1565 0.14379 0.1403
9LinNURBS 0.4713 0.1581 0.13883 0.135368 Fig. 12. Through-the-thickness distribution of interlaminar normal stress,
Closed Form 0.4678 0.1553 0.1430 0.1392 r zz ð0; 0; zÞ, in a simply-supported (SS2), antisymmetric (45/45) square laminated
plate, a/h = 50.
 H. Kapoor et al. / Composite Structures 106 (2013) 537–548 547

Table 8
Non-dimensionalized in-plane and interlaminar shear stresses in a simply-supported (SS1), sandwich square plate under sinusoidal loading.

a/h Source r xx r yy r xy r xz r yz GIS

(0) (0) (0.5)
4 9LinNURBS 0.8714 ) 0.15587 0.1157 – – 22
9LinNURBS 1.5896(0) 0.2651(0) 0.14397(0.2) – 33
9QuadNURBS 0.9167 0.16191 0.0939 0.26094 0.0682 22
9QuadNURBS 0.9522 0.1679 0.0977 0.2645 0.0683 33
16QuadNURBSKR – – – 0.2396(0.25) 0.1078(0.4) 33
16QuadNURBSKR – – – 0.2343(0.45) 0.1078(0.625) 44
16CubicNURBSKR – – – 0.2929 0.0781 33
16CubicNURBSKR – – – 0.2127 0.0789 44
9QuadFEM/Reddy 0.8918 0.1551 0.0901 0.2808 0.0746 33
3D Elasticity 1.556 0.2595 0.1481 0.239 0.1072

comparison of through thickness distribution of interlaminar nor- value of r xx ð0; 0; h=2Þ as 1.5896 as compared to 1.556 computed
mal stress with 3D elasticity and Kant’s solution. using 3D elasticity solution, 0.8918 computed using Reddy’s 4  4
mesh of 9QuadLagrange finite element and 0.9522 computed using
5.4. A simply-supported, anti-symmetric, angle-ply, square laminate 9QuadNURBS element with 2  2 mesh and 3  3 Gauss integration
scheme. Similar results are obtained for r  yy ð0; 0; h=2Þ as well. How-
A simply supported (SS2), 2 and 10 layered, (45–45), anti-sym- ever, in case of r  xy ða=2; b=2; h=2Þ, LinearNURBS element with knot
metric square laminate for various a/h ratios, subjected to sinusoi- insertion at {0.2} and with 3  3 Gauss integration scheme produces
dal loading is analyzed here. The material properties used here are 3D elasticity equivalent in-plane shear stress.
as follows: Interlaminar shear stresses are computed using 9QuadNURBS,
16QuadNURBSKR and 16CubicNURBS elements with 2  2 mesh
E1 =E2 ¼ 40; G12 =E2 ¼ 0:60; G23 =E2 ¼ 0:5; G13 ¼ G12 ;
and 3  3 and 4  4 Gauss integration schemes. k-refined NURBS
E2 ¼ E3 ; m12 ¼ m23 ¼ m13 ¼ 0:25 ð32Þ element produces interlaminar shear stresses which are closer to
Center deflection and in-plane stresses are computed using lin- 3D elasticity solution than 9QuadNURBS and 16CubicNURBS ele-
ear and quadratic NURBS element and are compared with closed- ments. Table 8 shows the comparison of in-plane and interlaminar
form solution. For this particular angle, in-plane stresses do not shear stresses with 3D elasticity and 4  4 mesh of 9QuadLagrange
vary with increasing a/h ratios. Table 7 shows comparison of center finite element.
deflection and in-plane stresses with closed form solution.
9QuadNURBS and 9LinNURBS elements with 2  2 mesh compares 6. Conclusions
well with the closed form solution.
Next, the interlaminar shear and normal stresses are computed This research paper describes the development of NURBS
for antisymmetric, two layered, (45/45), angle-ply laminate using Isogeometric finite element and post-processor for computation
quadratic and higher-order NURBS elements. In case of quadratic of interlaminar stresses in composite and sandwich plates. First-
NURBS element, 2  2 and 4  4 meshes are considered while for order, shear-deformable laminate composite plate theory is
higher-order NURBS elements, 2  2 mesh is considered. Fig. 9 utilized in deriving the governing equations using a variational
shows the through thickness distribution of r  xz for quadratic and formulation. Linear, quadratic, higher order and k-refined NURBS
higher order NURBS elements. Figs. 10–12 shows the through elements in parent domain are constructed. Lagrange finite
thickness distribution of r zz for higher order NURBS elements for element suffers from higher order stress gradient oscillations
various a/h ratios. Most of the elements converges to the same val- due to Gibbs phenomenon and require alternative stress recovery
ues of stresses. procedures. In this paper, a NURBS based post-processor is devel-
oped which computes interlaminar shear and normal stresses
5.5. A simply supported, sandwich plate directly from higher order gradients of NURBS basis in a single
step procedure. It is observed that linear NURBS element with
A simply supported (SS1) sandwich plate subjected to sinusoi- single knot insertion and higher order k-refined NURBS elements
dally distributed transverse loading is analyzed here. The face compute 3D elasticity equivalent in-plane and interlaminar stres-
sheet’s material (orthotropic) properties used here are as follows: ses respectively.

E1 ¼ 25E2 ; E2 ¼ 1E6; G12 ¼ 0:50E2 ; G23 ¼ 0:2;

Acknowledgments
G13 ¼ G12 ; m12 ¼ m13 ¼ 0:25 ð33Þ
Core material (transversely anisotropic) properties are as follows: I would like to thank ICTAS at Virginia Tech for funding this re-
search endevour. I would also like to thank Dr. Kapania and Dr.
E1 ¼ E2 ¼ 1E6; G13 ¼ G23 ¼ 0:06E6; G12 ¼ 0:16E6; Soni for their valuable discussions.
m12 ¼ 0:25 ð34Þ
References
Each face-sheet is assumed to have one tenth of the total thickness
of the sandwich plate. In-plane stresses are computed using [1] Piegl L, W Tiller. The NURBS book (Monograph in visual communication). 2nd
9QuadNURBS and 9LinNURBS elements, 2  2 mesh and different ed. Springer-Verleg; 1997.
[2] Farin G. Curves and surfaces. A practical guide. 5th ed. Morgan Kaufmann
Gauss integration schemes (GIS). Linear NURBS element with {0}
Publishers; 1993.
as single knot insertion computes in-plane stresses which are very [3] Cohen E, Lyche T, Reisenfield R. Discrete B-Splines and subdivision techniques
close to 3D elasticity solution. While 9QuadNURBS element behaves in computer aided geometric design and computer graphics. Comput Graph
Image Process 1980;14:87–111.
similar to 9QuadLagrange finite element with slightly better in-
[4] Rogers DF. Introduction to NURBS with historic prospective. Academic Press;
plane stresses. For example, 9LinNURBS element with {0} as single 1995.
knot insertion and with 3  3 Gauss integration scheme computes [5] Peters J, Reif U. Subdivision surface. Springer-Verlag; 1973.
548 H. Kapoor et al. / Composite Structures 106 (2013) 537–548

[6] Warren J, Weimer H. Introduction to NURBS with historic prospective: [31] Kapoor H, Kapania RK. Geometrically nonlinear NURBS Isogeometric finite
subdivision methods for geometric design. Morgan Kaufmann Publishers; element analysis of laminated composite plates. Compos Struct
1985. 2012;94(12):3434–47.
[7] Gordon W. Spline-blended interpolation through curve networks. J Meth Mech [32] Reddy JN. A simple higher order theory for laminated composite plates. J Appl
1969;18:931–52. Mech 1984;51:745–52.
[8] Gregory JA. N-sided surface patches. In: Gregory JA, editor. Mathematics of [33] Tessler A. Accurate interlaminar stress recovery from finite element analysis.
surfaces. Clarendon Press; 1965. NASA technical memorandum, NASA-TM-109149; September 1994.
[9] Loop CT, DeRose TD. A multisided generalization of Bezier surfaces. ACM Trans [34] Byun C, Kapania RK. Prediction of interlaminar stresses in laminated plates
Graph 1989;8:204–34. using the global orthogonal interpolation polynomials. AIAA J 1992;30(11).
[10] Bajaj C, Chen J, Xu G. Modeling with cubic a-patches. ACT Trans Graph [35] Lee K, Lee S. A post-processing approach to determine transverse stresses in
1995;14:103–33. geometrically nonlinear composite and sandwich structures. J Compos Mater
[11] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, 1996;37:2207–24.
NURBS, exact geometry and mesh refinement. Comput Meth Appl Mech Eng [36] Park JW, Kim YH. Displacement, stresses and their sensitivity coefficients in
2005;194:4135–95. composite panels. J Compos Mater 1999;33(13):1222–43.
[12] Kagan P, Fischer A, Yoseph PB. New B-Spline finite element approach for [37] Park BP, Park JW, Kim YH. Stress recovery in laminated composite and
geometrical design and mechanical analysis. Int J Numer Methods Eng sandwich panels undergoing finite rotation. Compos Struct 2003;59:227–35.
1998;41:435–58. [38] Noor AK, Kim YH, Peters JM. Transverse shear stresses and their sensitivity
[13] Kagan P, Fischer A, Yoseph PB. Mechanically based model: adaptive refinement coefficients in multi-layered composite panels. AIAA J 1994;32(6):1259–69.
of B-Spline finite element. Int J Numer Methods Eng 2003;57:1145–75. [39] Matsunaga H. Interlaminar stress analysis of laminated composite beams to
[14] Reissner E. On the theory of bending of elastic plates. J Math Phys global higher-order deformation theories. Compos Struct 2002;55:105–14.
1944;23:184–91. [40] Makeev A, Armanios EA. An iterative method for solving elasticity problems
[15] Reissner E. Effect of transverse shear deformation on the bending of elastic for composite laminates. ASME J Appl Mech 2000;67(1).
plates. J Appl Mech 1945;12:69–77. [41] Rolfes R, Rohwer K. Improved transverse shear stresses in composite finite
[16] Mindlin R. Influence of rotatory inertia and shear on flexural motions of elements based on the first order shear deformation theory. Int J Numer
isotropic, elastic plates. J Appl Mech 1951;18:31–8. Methods Eng 1997;40:51–60.
[17] Whitney J, Pagano N. Shear deformation in heterogeneous anisotropic plates. J [42] Kant T, Manjunatha BS. On accurate estimation of transverse stresses in
Appl Mech 1970;37:1031–6. multilayer composites. Comput Struct 1994;50(3):3365–451.
[18] Urthaler Y, Reddy JN. A mixed finite element for the nonlinear bending [43] Kant T, Gupta AB, Pendhari SS, Desai YM. Elasticity solution of cross-ply
analysis of laminated composite plates based on FSDT. Mech Adv Mater Struct composite and sandwich laminate. Comput Struct 2008;83(1):13–24.
2008;15:335–54. [44] Nosier A, Baharami A. Interlaminar stresses in antisymmetric angle-ply
[19] Chen S, Long YQ, Yao ZH, Chew SP. Application of the quadrilateral area co- laminates. Compos Struct 2007;78:18–33.
ordinate method: a new element for Mindlin–Reissner plate. Int J Numer [46] Barlow J. Optimal stress locations in finite element models. Int J Numer
Methods Eng 2006;66:1–45. Methods 1976;10(2):243–51.
[20] Kim KD, Chang SL, Sung CH. A 4-node corotational ANS shell element for [47] Oden JT, Brachli HJ. On the calculation of consistent stress distribution in finite
laminated composite structures. Compos Struct 2007;80(2):234–52. element approximation. Int J Numer Methods 1971;10(2):243–51.
[21] Mia-Duy N, Khennane A, Tran-Cong T. Computation of laminated composite [48] Hinton JS, Cambell E. Local and global smoothing of discontinuous finite
plates using integrated radial basis functions network. CMC: Comput, Mater element functions using a least square method. Int J Numer Methods
Continua 2007;5:63–77. 1974;8(3):461–80.
[22] Putcha NS, Reddy JN. A refined mixed shear flexible finite element for the [49] Chen DJ, Shah DK, Chan WS. Interfacial stress estimation using least-square
nonlinear analysis of laminated plates. Comput Struct 1986;22:529–38. extrapolation and local stress smoothing in a laminated composites. Comput
[23] Kant T, Kommineni JR. C0 finite element geometrically nonlinear analysis of Struct 1996;58:765–74.
fiber reinforced composite and sandwich laminates based on a higher-order [50] Zienkiewicz OC, Zhu JZ. The super convergent patch recovery and a posteriori
theory. Comput Struct 1992;45:511–20. error estimates. Part I: the recovery technique. Int J Numer Methods Eng
[24] Polit O, Touratier M. A multi layered/sandwich triangular finite element 1992;33(7):1331–64.
applied to linear and non-linear analyses. Compos Struct 2002:121–8. [52] Yuan J, Saleeb A, Gendy A. Projection stress: layerwise-equivalent formulation
[25] Phan ND, Reddy JN. Analysis of laminated composite plate using a higher order for accurate prediction of transverse stresses in laminated plate and shells. Int
shear deformation theory. Int J Numer Methods Eng 1985;21(12):2201–19. J Comput Eng Sci 1982;1(1):91–138.
[26] Ren JG, Hinton E. The finite element analysis of homogeneous and laminated [53] Prathap G, Naganarayana BP. Consistent thermal stress evaluation in finite
composite plates using a simple higher order theory. Commun Appl Numer elements. Comput Struct 1995;54(3):415–46.
Methods 1986;2(2):217–28.
[27] Reddy JN, Robbins Jr DH. Theories and computational models for composite
laminates. Appl Mech Rev 1994;47(6):147–69.
Further reading
[28] Carrera E. Historical review of zig–zag theories for multi layered plates and
shells. Appl Mech Rev 2003;56:287–303. [45] Senthil S, Batra RC. Interlaminar stresses in antisymmetric angle-ply
[29] Wang CM, Reddy JN. Shear deformable beams and plates: relationship with laminates. AIAA J 1999;37(11):1464–73.
classical solutions. Oxford (UK): Elsevier; 2000. [51] Yuan J, Saleeb A, Gendy A. Projection stress: layerwise-equivalent
[30] Reddy JN, Sandidge D. Mixed finite element models for laminated composite formulation for accurate prediction of transverse stresses in laminated
plates. J Eng Ind 1986;109:39–45. plate and shells. Int J Comput Eng Sci 1982;1(1):91–138.