# Analyzing Thermal Buckling in Curvilinearly Stiffened Composite Plates with
Arbitrary Shaped Cutouts Using Isogeometric Level Set Method

Balakrishnan Devarajan<EMAIL_ADDRESS>Department of Biomedical Engineering
and Mechanics, Virginia Polytechnic Institute and State University,
Blacksburg, VA 24061, USA

###### Abstract

In this paper we develop a new simple and effective isogeometric analysis for
modeling thermal buckling of stiffened laminated composite plates with cutouts
using level sets. We employ a first order shear deformation theory to
approximate the displacement field of the stiffeners and the plate. Numerical
modeling with a treatment of trimmed objects, such as internal cutouts in
terms of NURBS-based isogeometric analysis presents several challenges,
primarily due to need for using the tensor product of the NURBS basis
functions. Due to this feature, the refinement operations can only be
performed globally on the domain and not locally around the cutout. The new
approach can overcome the drawbacks in modeling complex geometries with
multiple-patches as the level sets are used to describe the internal cutouts;
while the numerical integration is used only inside the physical domain.
Results of parametric studies are presented which show the influence of ply
orientation, size and orientation of the cutout and the position and profile
of the curvilinear stiffeners. The numerical examples show high reliability
and efficiency of the present method compared with other published solutions
and ABAQUS.

###### keywords:

NURBS Isogeometric analysis, Laminated composite plate, First-order shear
deformation theory, Thermal buckling, Unitized structures, Stiffened composite
panel and Curvilinear composite stiffener

††journal: Journal of Composite Structures

##  1 Introduction

Due to a large requirement of laminates for a variety of engineering
applications, the use of plate structures with arbitrary cutouts are
inevitable. The presence of cutouts can significantly affect the behaviors of
structures. Studies on the analysis of eigenvalues and stability of laminates
with cutouts are of great importance for many practical applications including
airplane wing, fuselage and ribs.

In recent years, the isogeometric analysis (IGA) [1], [2] has become popular
since it inherently owns many great advantages including exact geometry
representation, higher-order continuity, simple mesh refinement, and avoiding
the traditional mesh generation procedure. The traditional mesh generation in
real world design problems could be complicated [3, 4, 5, 6, 4, 7, 8] Many
insights into mathematical properties [9], [10], integration method [11], [12]
and splines techniques [13] have been gained. The inherent characteristics of
IGA make it superior to the classical FEM in many aspects [1]. The IGA has
successfully applied to many engineering problems including plates and shells
[14], [15], [16], [17] fluid mechanics [9], damage and fracture mechanics
[18], contact mechanics [19], unsaturated flow problem in porous media [20],
and structural shape optimization [21].

The higher-order smoothness of the NURBS basis functions have huge benefits
with regards to the analysis of plates and shells. Based on the Kirchhoff
plate theory, construction of rotation-free isogeometric shell was first
proposed [22], and then developed for multiple NURBS patches by Kiendl [23].
Later, this element was extended for vibration and buckling analyses of
laminated plates [24], functionally graded plates [15], bending and buckling
analysis of laminated composite plates [17, 16, 25] and cloth simulation [19].

However, most of the problems investigated were simple single-patch geometry.
In modeling a more complex geometry like a plate with a clover shaped internal
cutout a trimmed NURBS surface can be useful. A NURBS surface and a set of
ordered boundary curves lying within the parameter space of the surface [26]
can be used to describe a trimmed NURBS surface. Shojaee et al. [24, 27]
divided a rectangular plate with a heart shaped cutout into multiple NURBS
patches, and to ensure continuity between adjacent patches, the authors
applied the bending strip method. However this approach is quite cumbersome to
implement and is also observed to complicate the unification of design and
analysis, which is the sole aim of IGA. Trimmed NURBS surfaces were handled by
Schmidt et al. [28] using a local reconstruction technique, which can handle
trimmed NURBS surfaces under an isogeometric analysis framework. The method
involved reconstructing trimmed NURBS patches using a least square
approximation or by an interpolation method having known the trimming curves.
The method is effective and can be easily applied to an arbitrarily shaped
area on the surface. However, the method could not exactly describe trimmed
boundaries and the number of sampling points heavily influenced the accuracy
of the reconstructed boundaries. In order to model complex geometry problems,
the IGA has been used in combination with an enrichment function similar to
XFEM, to develop XIGA [29]. It has also been used in combination with the
finite cell method (FCM) [30] and [31]. In this method, internal cutouts were
described using level sets. However, level set methods using IGA were
developed mostly for elastostatic problems, and are yet to be implemented for
thermal buckling analysis of stiffened laminated plates with complicated
internal cutouts. The objective of this research, is to fill this gap.

Thermal effects, especially thermal buckling is important to be considered in
design, analysis and weight optimization of aircraft structures[32]and [33].
Such problems have been researched by many in the past. For instance, Noor and
Burton [34],[35] and [36] used 3D elasticity to analyze the thermal buckling
of multilayered and sandwich plates. However, it is difficult to solve
problems which involve complicated geometry under different boundary
conditions. The problem becomes even more challenging if the plate is
reinforced with stiffeners.

Several authors in the past have studied the behavior of such stiffened
plates. A finite element algorithm was developed by Mukhopadhyay and Mukherjee
[37] to study the buckling of stiffened panels under uniaxial compression.
Prusty and Satsangi [38] used the same algorithm to perform structural
analysis of stiffened isotropic plates and shells. While the aforementioned
authors studied straight stiffeners, curvilinearly stiffened panels were
analyzed by Tamijani and Kapania [39], Zhao and Kapania [40], De and Kapania
[4, 41, 42, 43] and Shi et al. [44]. It has been found by Hao et al. [45] that
curvilinear stiffeners if used can improve the strength margins of panels
since they reinforce the cutouts in them. The tension field resulting due to
the presence of curvilinear stiffeners was examined and the loading path was
compared with results obtained using straight stiffeners.

In this work, thermal buckling analysis of stiffened laminated composite
plates with a cutout is presented. Level sets are used to describe the cutouts
and numerical integration is used only inside the physical domain. The plate
is assumed to be subjected to a uniform temperature rise for different
boundary conditions. The stiffness matrices and thermal force vector are
derived according to first-order shear deformation theory (FSDT). Firstly, the
present formulation is verified by comparison with available references from
literature. The significant effects of aspect ratio, length-to-thickness
ratio, fiber orientation and boundary condition on the behavior of composite
plates are then analyzed and useful conclusions are derived from the numerical
results.

##  2 Mathematical foundations

###  2.1 Laminated composite panel

According to the first-order shear deformation theory (FSDT) the displacement
field is considered as the first-order Taylor series expansion of the mid-
plane displacement variables with respect to plate thickness as follows:
Consider a stiffened composite plate with a circular cutout as shown in Figure
1. The mid-plane of the panel $O_{xy}$ is used as the reference plane for the
global coordinate system. The plate is assumed to have a uniform thickness
with no plydrops. Let $u$ be the displacements in the $x$-axis and $v$ be the
in-plane displacements along the $y$-axis. Let $w$ be the transverse
deflection along the $z$-axis.

Figure 1: Geometry and nomenclature of a stiffened composite panel with a
central cutout.

Let $\beta_{y}$ and $\beta_{x}$ be rotation components of the panel around the
$x$\- and $y$-axes, respectively. The displacement components of the panel can
now be written as,

$\begin{array}[]{cc}&u\left(x\text{,}y\text{,}z\text{,}t\right)=u_{0}\left(x\text{,}y\text{,}z\text{,}t\right)+z\beta_{x}\left(x\text{,}y\text{,}t\right)\\\
&v\left(x\text{,}y\text{,}z\text{,}t\right)=v_{0}\left(x\text{,}y\text{,}z\text{,}t\right)+z\beta_{y}\left(x\text{,}y\text{,}t\right)\\\
&w\left(x\text{,}y\text{,}z\text{,}t\right)=w_{0}\left(x\text{,}y\text{,}z\text{,}t\right)\end{array}$
(1)

The plate strain energy $U_{p}$ can be written as,

$U_{p}=\frac{1}{2}\underset{\Omega}{\hskip 2.79999pt\int\int\hskip
2.79999pt}\boldsymbol{\varepsilon}_{\mathbf{p}}^{\mathbf{L}^{\mathrm{T}}}{\mathbf{D}}_{\mathbf{p}}\boldsymbol{\varepsilon}_{\mathbf{p}}^{\mathbf{L}}\mathrm{d}\Omega$
(2)

Where, $\mathbf{D}_{p}$ depends on the material property, stacking sequence,
ply orientation and thickness of the plate. Derivation of $\mathbf{D}_{p}$ is
explained in Section 2.2. The generalized strains
$\boldsymbol{\varepsilon}_{\mathbf{p}}^{\mathbf{L}}$ can be expressed in terms
of the generalized displacements $\mathbf{u_{p}}$ of the panel as,

$\boldsymbol{\varepsilon}_{\mathbf{p}}^{\mathbf{L}}=\begin{Bmatrix}\varepsilon_{x}^{0}\\\
\varepsilon_{y}^{0}\\\ \gamma_{xy}^{0}\\\ \kappa_{x}^{0}\\\ \kappa_{y}^{0}\\\
\kappa_{xy}^{0}\\\ \gamma_{xz}^{0}\\\ \gamma_{yz}^{0}\\\
\end{Bmatrix}=\begin{bmatrix}\frac{\partial}{\partial{x}}&0&0&0&0\\\
0&\frac{\partial}{\partial{x}}&0&0&0\\\
\frac{\partial}{\partial{x}}&\frac{\partial}{\partial{x}}&0&0&0\\\
0&0&0&\frac{\partial}{\partial{x}}&0\\\
0&0&0&0&\frac{\partial}{\partial{x}}\\\
0&0&0&\frac{\partial}{\partial{x}}&\frac{\partial}{\partial{x}}\\\
0&0&\frac{\partial}{\partial{x}}&1&0\\\
0&0&\frac{\partial}{\partial{x}}&0&1\\\ \end{bmatrix}\begin{Bmatrix}u_{0}\\\
v_{0}\\\ w_{0}\\\ \beta_{x}\\\
\beta_{y}\end{Bmatrix}=\mathbf{B_{p}}\mathbf{u_{p}}$ (3)

The strain energy of the composite panel $U_{p}$ is,

$U_{p}=\frac{1}{2}\iint_{\Omega}\mathbf{u_{p}}^{T}\mathbf{B_{p}}^{T}\mathbf{D_{p}}\mathbf{B_{p}}\mathbf{u_{p}}d\Omega$
(4)

For thermal bucking analysis, the geometric stiffness due to inplane stress
resultants is considered. Thus, the potential $W_{p}$ of the panel due to in-
plane stress vector
${\overline{\sigma}}_{\mathbf{p}}={\left(\sigma_{x}^{0}\hskip 2.79999pt\hskip
2.79999pt\hskip 2.79999pt\sigma_{y}^{0}\hskip 2.79999pt\hskip 2.79999pt\hskip
2.79999pt\tau_{xy}^{0}\hskip 2.79999pt\hskip 2.79999pt\hskip 2.79999pt0\hskip
2.79999pt\hskip 2.79999pt\hskip 2.79999pt0\right)}$ is,

$W_{p}=\frac{1}{2}\underset{\Omega}{\hskip 2.79999pt\int\int\hskip
2.79999pt}{\overline{\sigma}}_{\mathbf{p}}\boldsymbol{\varepsilon}_{\mathbf{p}}^{\mathbf{NL}}\mathrm{d}\Omega$
(5)

The nonlinear strain factor matrix is:

$\left[\mathbf{B}_{p}^{\mathit{NL}}\right]=\begin{bmatrix}0&0&\frac{\partial}{\partial
x}&0&0\\\ 0&0&\frac{\partial}{\partial y}&0&0\\\
0&0&0&\frac{\partial}{\partial x}&0\\\ 0&0&0&\frac{\partial}{\partial y}&0\\\
0&0&0&0&\frac{\partial}{\partial x}\\\ 0&0&0&0&\frac{\partial}{\partial
y}\end{bmatrix}$ (6)

The in-plane stress matrix can be written as:

${\boldsymbol{\sigma}}_{\mathbf{p}}={\left(\begin{array}[]{cccccc}t\sigma_{y}^{0}&t\tau_{xy}^{0}&0&0&0&0\\\
t\tau_{xy}^{0}&t\sigma_{y}^{0}&0&0&0&0\\\
0&0&\frac{t^{3}}{12}\sigma_{y}^{0}&\frac{t^{3}}{12}\tau_{xy}^{0}&0&0\\\
0&0&\frac{t^{3}}{12}\tau_{xy}^{0}&\hskip
2.79999pt\frac{t^{3}}{12}\sigma_{y}^{0}&0&0\\\
0&0&0&0&\frac{t^{3}}{12}\sigma_{y}^{0}&\frac{t^{3}}{12}\tau_{xy}^{0}\\\
0&0&0&0&\frac{t^{3}}{12}\tau_{xy}^{0}&\hskip
2.79999pt\frac{t^{3}}{12}\sigma_{y}^{0}\end{array}\right)},$ (7)

Hence, Equation 5 can be written as

$W_{p}=\frac{1}{2}\underset{\Omega}{\hskip 2.79999pt\int\int\hskip
2.79999pt}\mathbf{u}_{\mathbf{p}}^{\mathbf{T}}\mathbf{B}_{\mathbf{p}}^{\mathbf{NL}^{\mathrm{T}}}{\mathbf{\sigma}}_{\mathbf{p}}\mathbf{B}_{\mathbf{p}}^{\mathbf{NL}}{\mathbf{u}}_{\mathbf{p}}\mathrm{d}\Omega$
(8)

###  2.2 Orthotropic layers

The constitutive equation in the local coordinate system for the $k^{th}$
orthotropic elastic lamina can be derived from Hooke's law as

$\begin{Bmatrix}\sigma_{1}\\\ \sigma_{2}\\\ \tau_{12}\\\ \sigma_{3}\\\
\tau_{31}\\\
\tau_{23}\end{Bmatrix}^{(k)}=\begin{bmatrix}Q_{11}&Q_{12}&0&0&0&0\\\
Q_{21}&Q_{22}&0&0&0&0\\\ 0&0&Q_{66}&0&0&0\\\ Q_{31}&Q_{32}&0&Q_{33}&0&0\\\
0&0&0&0&Q_{55}&0\\\
0&0&0&0&0&Q_{44}\end{bmatrix}^{(k)}\begin{Bmatrix}\varepsilon_{1}-\alpha_{1}\mathrm{\Delta}T\\\
\varepsilon_{2}-\alpha_{2}\mathrm{\Delta}T\\\ \gamma_{12}\\\
\varepsilon_{3}-\alpha_{3}\mathrm{\Delta}T\\\ \gamma_{31}\\\
\gamma_{23}\end{Bmatrix}^{(k)}$ (9)

where $Q_{ij}$ are the elastic coefficients in the material coordinate system,
$T$ is the temperature change and $\alpha_{i}$ is the thermal coefficient of
expansion in the principal $i^{th}$-direction. The transformed material
constants are expressed as,

$\begin{Bmatrix}{\overline{Q}}_{11}\\\ {\overline{Q}}_{12}\\\
{\overline{Q}}_{22}\\\ {\overline{Q}}_{16}\\\ {\overline{Q}}_{26}\\\
{\overline{Q}}_{66}\end{Bmatrix}=\begin{bmatrix}c^{4}&2c^{2}s^{2}&s^{4}&4c^{2}s^{2}\\\
c^{2}s^{2}&c^{4}+s^{4}&c^{2}s^{2}&-4c^{2}s^{2}\\\
s^{4}&2c^{2}s^{2}&c^{4}&4c^{2}s^{2}\\\
c^{3}s&\mathit{cs}(s^{2}-c^{2})&-\mathit{cs}^{3}&2\mathit{cs}(s^{2}-c^{2})\\\
\mathit{cs}^{3}&\mathit{cs}(c^{2}-s^{2})&-c^{3}s&2\mathit{cs}(c^{2}-s^{2})\\\
c^{2}s^{2}&-2c^{2}s^{2}&c^{2}s^{2}&{(c^{2}-s^{2})}^{2}\end{bmatrix}\begin{Bmatrix}Q_{11}\\\
Q_{12}\\\ Q_{22}\\\ Q_{66}\end{Bmatrix}\text{,}$ (10)
$\begin{bmatrix}{\overline{Q}}_{44}&{\overline{Q}}_{13}\\\
{\overline{Q}}_{45}&-{\overline{Q}}_{36}\\\
{\overline{Q}}_{55}&{\overline{Q}}_{23}\end{bmatrix}=\begin{bmatrix}c^{2}&s^{2}\\\
-\mathit{cs}&\mathit{cs}\\\
s^{2}&c^{2}\end{bmatrix}\begin{bmatrix}Q_{44}&Q_{13}\\\
Q_{55}&Q_{23}\end{bmatrix}\text{,}\hskip 3.50006pt{\overline{Q}}_{33}=Q_{33}$
(11) $\begin{Bmatrix}\alpha_{x}\\\ \alpha_{y}\\\
\alpha_{\mathit{xy}}\end{Bmatrix}=\begin{bmatrix}c^{2}&s^{2}\\\ s^{2}&c^{2}\\\
2\mathit{cs}&-2\mathit{cs}\end{bmatrix}\begin{Bmatrix}\alpha_{1}\\\
\alpha_{2}\end{Bmatrix}\text{,}\hskip 3.50006pt\alpha_{z}=\alpha_{3}$ (12)

where $c=\cos\theta$ and $s=\sin\theta$ .
Using this, $\mathbf{D}_{p}$ can now be written as

$\mathbf{D}_{p}=\begin{bmatrix}\mathbf{A}&\mathbf{B}&\mathbf{0}\\\
\mathbf{B}&\mathbf{D}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{A}^{\text{s}}\end{bmatrix}$ (13)

where,

$(A_{\mathit{ij}}\text{,}{B}_{\mathit{ij}}\text{,}D_{\mathit{ij}})=\int_{-h/2}^{h/2}\overline{Q}_{\mathit{ij}}(1\text{,}z\text{,}z^{2})\text{d}z\text{,}\hskip
10.00002ptA_{\mathit{ij}}^{\text{s}}=K\int_{-h/2}^{h/2}\overline{Q}_{\mathit{ij}}\hskip
1.19995pt\text{d}z$ (14)

where $A_{\mathit{ij}}$, $B_{\mathit{ij}}$ and $D_{\mathit{ij}}$ are valid for
$i,j=1,2,6$, and $A_{\mathit{ij}}^{\text{s}}$ for $i,j=4$, 5 according to the
Voigt notation. $K$ denotes the transverse shear correction coefficient. A
value of $K=5/6$ was used for the analyses.

###  2.3 Curvilinear stiffener

Consider a stiffener attached to a plate as shown in Figure 2.

Figure 2: Composite plate (red) stiffened by a rectangular stiffener (gray)

The width and height of the rectangular stiffener are given as $b_{s}$ and
$h_{s}$, respectively. Details of the coordinate system and nomenclature can
be found in [40]. The stiffener is modeled using 1D $2^{nd}$ order NURBS basis
function. The stiffener is modeled as a Timoshenko beam. A detailed
description of how displacement compatibility between the control points of
the stiffener and the plate was achieved can be found in [46] and [47].

###  2.4 Isogeometric analysis

The concept of IGA is to use the same basis functions to generate CAD and to
discretize the solution to the boundary value problem. This isoparametric
formulation in combination with the properties of NURBS basis functions allows
the geometry to be modeled exactly from the coarsest level of refinement, and
maintains this exact representation as the mesh is refined. The elements are
constructed by a knot span which is an array of non-decreasing sequence of
parameter values [1]. The NURBS basis functions are used to construct the
exact geometry, as well as the corresponding solution space. In this section,
we present a primer on NURBS and also some basic formulations of IGA. Knot
vectors can be defined as:

$\Xi=\\{\xi_{1},\xi_{2},...,\xi_{n+p+1}\\}$ (15)

B-spline basis functions are defined recursively starting with piecewise
constants $p=0$

$N_{i}^{0}(\xi)=\left\\{\begin{array}[]{ccc}1&\mathit{if}&\xi_{i}<\xi<\xi_{i+1}\\\
0&&\text{otherwise}\end{array}\right.$ (16)

For higher values like $p=0,1,2...$ the basis functions are defined using a
recursion relation

$N_{i}^{p}(\xi)=\frac{\xi-\xi_{i}}{\xi_{i+p}-\xi_{i}}N_{i}^{p-1}(\xi)+\frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}}N_{i+1}^{p-1}(\xi)$
(17)

where $\xi_{1}$ is the $i^{th}$ knot, $n$ is the number of basis functions and
$p$ is the polynomial order.
We denote Non-uniform Rational B-splines (NURBS) functions expressed as

$R_{i,p}(\xi)=\frac{N_{i,p}(\xi)w_{i}}{\sum_{j=1}^{n}N_{j,p}(\xi)w_{j}}$ (18)

Where $w_{i}$ is the weight corresponding to the $i^{th}$ control point. For a
2-dimensional surface, the basis functions are expressed as:

$R_{i\text{,}j}^{p\text{,}q}(\xi\text{,}\eta)=\frac{N_{i\text{,}p}(\xi)N_{j\text{,}q}(\eta)w_{i\text{,}j}}{\sum_{i=1}^{n}\sum_{j=1}^{m}N_{i\text{,}p}(\xi)N_{j\text{,}q}(\eta)w_{i\text{,}j}}=\frac{N_{i\text{,}p}(\xi)N_{j\text{,}q}(\eta)w_{i\text{,}j}}{W(\xi\text{,}\eta)}$
(19)

The 2D surface is generated using the basis functions, $R_{i,j}$, and the
control points $P_{i,j}$

$\boldsymbol{S}(\xi\text{,}\eta)=\sum_{i=1}^{n}\sum_{j=1}^{m}R_{i\text{,}j}^{p\text{,}q}(\xi\text{,}\eta){\boldsymbol{P}}_{i\text{,}j}\text{,}$
(20)

Further details can be found in [1] and [25].

###  2.5 Level Set Function

A signed distance function is utilized here to identify the enriched elements
(elements cut by the cutout boundary), the outer elements (elements that are
part of the plate) and inner elements (elements that fall completely within
the cutout geometry) using the equation of a circle, given as:

$\phi(\mathbf{X})=|\mathbf{X}-\mathbf{X_{c}}|-r$ (21)

Where $\phi(\mathbf{X})$ is the level set value at coordinates $\mathbf{X}$.
$\mathbf{X_{c}}$ is the coordinates of the center of the circular cutout and
$r$ is the radius of the cutout. These values are stored at all control points
for later use. Next, the coordinates of the points where the geometry of the
cutout intersects the element is computed (Figure 3). This process will
approximate the shape of the cutout as a polygon with sides equal to the
number of enriched elements. This process has no impact on elements that are
not cut by the interface.

The implemented enrichment function is:

$\psi(\mathbf{X})=\sum_{I=1}^{4}|\phi_{I}|N(\mathbf{X})_{I}-|\sum_{I=1}^{4}\phi_{I}N(\mathbf{X})_{I}|$
(22)

The enrichment function $\psi$ at coordinates $\mathbf{X}$ is evaluated as a
combination of control point level set values and the shape functions.
Compared to other choices for enrichment functions this one has two major
advantages. As the enrichment function has a zero value outside the cut
elements as well as in the cut element control points there will be no
contribution to the global stiffness matrix from an element that are not cut
by the interface and the displacements evaluated at the control points will
not have any extra contribution either.

###  2.6 Triangulation and Numerical Integration

Figure 3 shows a non enriched element (A) and an enriched element (B) in the
$x-y$ coordinate system. Since some part of such enriched elements lie within
the cutout, for numerical integration purposes, they have to be treated
differently. For this, the coordinates of the points where the geometry of the
cutout intersects the element is computed (Figure 3). The creation of
additional points will enable the division of such a quadrilateral element
into triangles and the process will approximate the shape of the cutout as a
polygon with sides equal to the number of enriched elements

Consider one such triangular geometry ($\Delta_{1}$). Let the function
$f(x,y)$ determine some property at $[x,y]$, a point within this geometry.
Suppose another property, that is evaluated by the integration of $f(x,y)$
over this geometry is of interest. An analytical calculation of the integrand
is either impossible or very expensive. It is convenient to use Gauss
quadrature formula that approximate the integral as a weighted sum of the
function values at a number of predetermined points. Note that the Gauss
integration computed using $n$ points is exact for a function of order $2n-1$
in 1D. These points are given in a transformed coordinate system (natural
coordinate system).
Thus, to integrate $f(x,y)$ over $\Delta_{1}$,

$I=\int_{\Delta_{1}}f(x,y)dxdy$ (23)

a transformation into this natural coordinate system and a change of variables
is necessary. This transformation is performed using a Jacobian matrix,
$\mathbf{J_{el}}$ and is defined by

$\begin{bmatrix}dx\\\ dy\end{bmatrix}=\mathbf{J_{el}}\begin{bmatrix}d\xi\\\
d\eta\end{bmatrix}$ (24)

Where,

$\mathbf{J_{el}}=\begin{bmatrix}\frac{\partial x}{\partial\xi}&\frac{\partial
x}{\partial\eta}\\\ \frac{\partial y}{\partial\xi}&\frac{\partial
y}{\partial\eta}\end{bmatrix}$ (25)

Then it follows that

$I=\frac{1}{2}\int_{0}^{1}\int_{0}^{1-\eta}f(x(\xi,\eta),y(\xi,\eta))\begin{vmatrix}\frac{\partial
x}{\partial\xi}&\frac{\partial x}{\partial\eta}\\\ \frac{\partial
y}{\partial\xi}&\frac{\partial y}{\partial\eta}\end{vmatrix}d\xi d\eta$ (26)

By numerical integration the quantity can be evaluated

$I\cong\frac{1}{2}\sum_{1}^{i}w_{i}f(x(\xi,\eta),y(\xi,\eta))\begin{vmatrix}\frac{\partial
x}{\partial\xi}&\frac{\partial x}{\partial\eta}\\\ \frac{\partial
y}{\partial\xi}&\frac{\partial y}{\partial\eta}\end{vmatrix}$ (27)

Where $i$ is the number of Gauss points in the triangle. This will be useful
in codes where element stiffness matrices for elements with enriched nodes are
evaluated.

Figure 3:  Mappings and transformation from physical space to natural space

The element stiffness matrix will have a shape as follows:

$\mathbf{K_{e}^{ij}}=\iint_{elem}\mathbf{B_{i}^{T}}\mathbf{D}\mathbf{B_{j}}tdxdy$
(28)

Where,

$\mathbf{B_{i}}=[\mathbf{B_{STD}}\ |\ \mathbf{B_{ENR}}]$ (29)

Where $\mathbf{B_{i}}$ could be one of the three matrices, $\mathbf{B_{p}}$,
$\mathbf{B_{s}}$ or $\mathbf{B_{p}^{NL}}$. The subscripts $\mathbf{STD}$ and
$\mathbf{ENR}$ refer to standard non-enriched elements and enriched elements
respectively. For a four noded isoparametric element the shape functions are
defined in the isoparametric domain. However, the $\mathbf{B}$ matrices are
functions of the shape functions which are defined in the isoparametric
coordinate system $\xi-\eta$ (Figure 4) , but the Gauss points and weights
associated for a triangular domain aren’t.

Figure 4: Element and triangle transformed into coordinate system appropriate
for $B_{matrix}$ calculation

They are defined in the $\xi^{\prime}-\eta^{\prime}$ coordinate system (Figure
5 with the vertices of the triangle located at:

$\left[\xi{}_{1}^{\prime},\eta_{1}^{\prime}\right]=\left[0\
0\right];\left[\xi{}_{2}^{\prime},\eta_{2}^{\prime}\right]=\left[1\
0\right];\left[\xi{}_{1}^{\prime},\eta_{1}^{\prime}\right]=\left[0\ 1\right];$
(30)

Hence, the mapping of these Gauss points to the isoparametric coordinate
system is to be performed for each triangle.

These Gauss points can be mapped using nodal shape functions and isoparametric
coordinates of vertices as follows:

$\begin{bmatrix}\xi_{gp}&\eta_{gp}\end{bmatrix}_{\Delta_{1}}=\begin{bmatrix}1-\xi{}^{{}^{\prime}}_{gp}-\eta{}^{{}^{\prime}}_{gp}\\\
\xi{}^{{}^{\prime}}_{gp}\\\
\eta{}^{{}^{\prime}}_{gp}\end{bmatrix}^{T}_{\Delta_{1}}\begin{bmatrix}\xi_{1}&\eta_{1}\\\
\xi_{2}&\eta_{2}\\\ \xi_{3}&\eta_{3}\\\ \end{bmatrix}_{\Delta_{1}}$ (31)

The weight associated with each Gauss point will be scaled using the ratio of
the areas of the triangles in the $\xi$-$\eta$ coordinate system and the
$\xi^{\prime}$-$\eta^{\prime}$ coordinate system. We now have Gauss points and
associated weights in the isoparametric coordinate system to in order to
perform numerical integration.

(a) 3 Gauss points

(b) 7 Gauss points

Figure 5: Triangle transformed into coordinate system compatible with given
Gauss-points

Now, all the tools needed for the numerical integration are available. The
stiffness matrix contribution from the triangle can be calculated now as

$\iint_{\Delta}\mathbf{B_{i}^{T}}(x,y)\mathbf{D}\mathbf{B_{j}}(x,y)tdxdy=\frac{1}{2}\sum_{k=1}^{n}\mathbf{B_{i}}(x,y)^{T}\mathbf{D}\mathbf{B}_{j}(x,y)tw_{k}\left|\mathbf{J}_{\Delta}\right|\left|\mathbf{J}_{el}\right|$
(32)

Where $B_{i}$ and $B_{j}$ are either $\mathbf{B_{STD}}$ or $\mathbf{B_{ENR}}$,
depending on what part of the stiffness matrix are currently being calculated
and $n$ is the number of Gauss points in the element.

The target expression (a sample entry in a $\mathbf{B_{STD}}$ matrix) is,

$N_{1,x}=\frac{\partial\xi}{\partial x}\frac{\partial
N_{1}}{\partial\xi}+\frac{\partial\eta}{\partial x}\frac{\partial
N_{1}}{\partial\eta}$ (33)

Where $N$ is the basis function. While, for the enriched B-matrix
$\mathbf{B_{ENR}}$,

$\frac{\partial(\psi N_{1})}{\partial x}=N_{1}\frac{\partial\psi}{\partial
x}+\frac{\partial N_{1}}{\partial x}\psi$ (34)

Where, $\psi$ is the enrichment function.

###  2.7 A Three point parametrization of the curvilinear stiffener

Previous researchers [39], [40] and [44] modeled curvilinear stiffeners using
a parabolic curve. Using the parametric equation of the parabola and the start
and the end point coordinates, one can obtain the three control points which
could define the Bezier curve [48]. Since any parabolic stiffener can be
represented using three control points, the coordinates of the end control
points $\Delta\epsilon$ and Point 2 (where $\delta_{dist}$ refers to the point
[$\delta_{dist}$ ,$\delta_{dist}$]) were used to represent the curvilinear
stiffener.

Figure 6: Parametrization of curvilinear stiffener

##  3 Results and Discussion

To demonstrate and validate the method, thermal buckling analysis of
curvilinearly stiffened laminated composite plates with various cutouts is
carried out. Results are compared with existing literature for some cases
while for some others, using commercially popular software like ABAQUS.
Parametric studies were performed which show the influence of both the cutout
profile and the shape and position of the curvilinear stiffeners in the
results.

###  3.1 Thermal buckling of laminated composite plates

Over the years, and most recently [49], several researchers have analyzed the
thermal buckling of symmetric laminated plates. A four-layer
${[}0^{\circ}/90^{\circ}/90^{\circ}/0^{\circ}{]}$ laminated plate with the two
ratios of $a/h=10$ (Figure 2) and $100$ were investigated for [50] two sets of
boundary conditions. The plate is assumed to be a square of edge length 1 unit
to make normalizations easy.

$CCCC:\hskip 3.50006ptu_{0}=0,v_{0}=0,\psi_{x}=0,\psi_{y}=0,w_{0}=0$ (35)

on all four edges.

$\mathit{SSSS}:\left\\{\begin{array}[]{c}u_{0}=v_{0}=w_{0}=\phi_{y}=\phi_{z}=0\hskip
10.00002pt\text{on}\hskip 3.50006ptx\hskip 3.50006pt=\hskip
3.50006pt0\text{,}\hskip 3.50006pta\\\
u_{0}=v_{0}=w_{0}=\phi_{x}=\phi_{z}=0\hskip 10.00002pt\text{on}\hskip
3.50006pty\hskip 3.50006pt=\hskip 3.50006pt0\text{,}\hskip
3.50006ptb\end{array}\right.$ (36)

Material properties are defined as [51], [50] and[49]

$\begin{array}[]{c}E_{L}/E_{T}=15\text{,}\hskip
3.50006ptG_{\mathit{LT}}/E_{T}=0.5\text{,}\hskip
3.50006ptG_{\mathit{TT}}/E_{T}=0.3356\text{,}\\\
\nu_{\mathit{LT}}=0.3\text{,}\hskip
3.50006pt\nu_{\mathit{TT}}=0.49\text{,}\hskip
3.50006pt\alpha_{L}/\alpha_{0}=0.015\text{,}\hskip
3.50006pt\alpha_{T}/\alpha_{0}=1\end{array}$ (37)

where $L$ and $T$ refer to the directions parallel and perpendicular to the
fibers and $\alpha_{0}$ is the normalization factor for the thermal expansion
coefficient. Table LABEL:tab:ls_compositenostiff illustrates the comparison of
critical temperature obtained by NURBS-based IGA and FEM-Q16 [51], [50]
and[49]. The critical buckling temperature values of the thinner plate have
been scaled by 100.

###  3.2 Isotropic plate with a circular hole at the center

In order to illustrate the performance of IGA code while modeling plates with
holes, an isotropic plate with a circular hole at the center was considered.
Variation of the critical thermal load with respect to the radius of the hole
was compared against Avci et al.[52]. The results were obtained for simply
supported (SSSS) and clamped (CCCC) boundary conditions. The material
properties of the plate are given as:

$\begin{array}[]{c}E\hskip 3.50006pt=\hskip 3.50006pt208\hskip
3.50006ptGPa\text{,}\hskip 3.50006pt\nu=\hskip 3.50006pt0.3\text{,}\hskip
3.50006pt\alpha=\hskip 3.50006pt1.17e-5\hskip
3.50006pt^{\circ}C^{-1}\end{array}$ (38)

As seen from the Table 1, the results are in good agreement with [52].

Table 1: Variation of critical load with respect to radius of the circular cutout | CCCC | SSSS
---|---|---
$\frac{diameter}{width}$ | Present | Avci | Present | Avci
0 | 29.35 | 29.16 | 10.83 | 10.97
0.1 | 28.01 | 28.57 | 10.64 | 10.82
0.2 | 29.17 | 29.07 | 10.52 | 10.68
0.3 | 35.27 | 35.37 | 11.01 | 11.14
0.4 | 49.53 | 49.67 | 12.48 | 12.57
0.5 | 78.74 | 78.55 | 15.41 | 15.34

###  3.3 Composite plate with elliptical cutout

To demonstrate the versatility and to validate the method, buckling analysis
of plates with an elliptical cutout was considered (Figure 7). Results are
compared with those obtained using ABAQUS. A unit square domain is considered
with an elliptical cutout. The selected examples are for comparison purpose
only and do not represent the full capabilities of the proposed method.

Figure 7: Dimensions and orientation of an elliptical cutout

The ellipse is defined by the semi major axis, the semi minor axis and the
angle of orientation $\theta$. For validation purposes, an ellipse of semi-
major axis = 0.2, semi-minor axis = 0.1 and $\theta=0$ was considered.

Table 2: The convergence of the critical thermal buckling load of a four-layer $[0^{\circ}/90^{\circ}/90^{\circ}/0^{\circ}]$ laminated composite square plate. Number of elements | ABAQUS | Present
---|---|---
16 | 0.636 | 0.653
64 | 0.406 | 0.421
256 | 0.3782 | 0.401
1024 | 0.377 | 0.381

It can be observed from Table 2 that level set IGA approach yields relatively
the same rate of convergence as ABAQUS. For mode shape comparison (Figure 8)
purposes, a plate with an elliptical cutout ($\theta=45^{\circ}$) was modeled
in ABAQUS with an irregular mesh composed of 1028 elements. S8R: An 8-node
doubly curved thick shell, reduced integration element was used to mesh the
domain. IGA mode shape plot was obtained using 1024 elements.

(a) ABAQUS

(b) IGA

Figure 8: Buckling mode shape plots of a plate with an elliptical cutout
($\theta=45^{\circ}$)

Next, the variation of critical buckling load with respect to change in
eccentricity and angle of orientation of the cutout was analyzed. Table 3
presents the results for three different eccentricities. The results of the
present method (LSM) are compared with those obtained using the Multi-Patch
method (MPM)

Table 3: Variation of thermal buckling load of a four-layer $[0^{\circ}/90^{\circ}/90^{\circ}/0^{\circ}]$ laminated composite square plate. $\theta$ |  | Semi-minor axis=0.1
---
Semi-major axis=0.2
| Semi-minor axis=0.1
---
Semi-major axis=0.3
| Semi-minor axis=0.2
---
Semi-major axis=0.3
| LSM | MPM | LSM | MPM | LSM | MPM
0 | 0.381 | 0.379 | 0.413 | 0.402 | 0.633 | 0.633
15 | 0.395 | 0.390 | 0.443 | 0.438 | 0.659 | 0.653
30 | 0.427 | 0.421 | 0.532 | 0.525 | 0.729 | 0.715
45 | 0.469 | 0.460 | 0.641 | 0.630 | 0.809 | 0.807

###  3.4 Curvilinearly stiffened composite panels with a noncentric circular
cutout

In this sub section, the influence of the noncentric circular cutouts on
clamped stiffened composite panels is studied. The stiffener is isotropic with
Young’s modulus = $E_{T}$, Poisson’s ratio = $\nu_{LT}$ and coefficient of
thermal expansion = $\alpha_{0}$. The stiffness ratio
$\gamma=\mathit{EI}/\mathit{bD}$ and the area ratio
$\delta=A_{s}/{\mathit{bt}}_{p}$ are 5 and 0.1 respectively unless specified.
To check the robustness of the level set method, results for clamped
curvilinearly stiffened composite plate are compared with those obtained using
ABAQUS. A noncentric circular cutout of radius $0.15$ with center at
($0.3,0.7$) was assumed. The curvilinear stiffener configuration was adopted
from [39] and scaled proportionally to match the dimensions of the plate. From
the ABAQUS library, the plate was meshed using 3040 S8R elements and
curvilinear stiffener was meshed using 192 elements S8R elements (48 divisions
along the arc length and 4 divisions along the depth direction). The IGA plate
model was meshed using 4096 elements and the curvilinear stiffener using 32
elements. The critical buckling temperature result was observed to converge
for this mesh configuration. The first five eigen buckling models were
extracted for each case. As seen from Figure 9 the buckling mode shapes
obtained using IGA and ABAQUS are comparable.

| Present | ABAQUS
---|---|---
|  |  |  |
Mode |  |  $\Delta T_{cr}$ |  |  $\Delta T_{cr}$
1 |  |  0.426 |  |  0.420
2 |  |  0.535 |  |  0.523
3 |  |  0.708 |  |  0.698
4 |  |  0.910 |  |  0.906
5 |  |  1.070 |  |  1.007
Figure 9: First five eigenmode shape plots for a plate with curvilinear
stiffener using LSM and ABAQUS, a commercial available software.

###  3.5 Curvilinearly stiffened composite panels with an elliptical central
cutout

In this section, results of curvilinearly stiffened composite panels with
elliptical cutouts are presented. An ellipse of semi-major axis = 0.2 and
semi-minor axis = 0.1 is considered for all cases. The variation of the
critical thermal load with respect to angle orientations, ply orientations,
stiffness ratios and different stiffener profiles (defined by $\Delta\epsilon$
Figure 6) are presented in Table 5 and Table 6. The results obtained using
present method (LSM) and the Multi-Patch method (MPM) are in excellent
agreement.

Table 5: Variation of critical buckling load with respect to angle of
orientation of the elliptical cutout for different ply orientations for
$\gamma=5$ $\Delta\epsilon=0$
---
$\theta$ |  | Antisymmetric
---
Cross-Ply
| Symmetric
---
Cross-Ply
| Symmetric
---
Angle-Ply
| Antisymmetric
---
Angle-Ply
| LSM | MPM | LSM | MPM | LSM | MPM | LSM | MPM
0 | 0.475 | 0.473 | 0.413 | 0.406 | 0.425 | 0.423 | 0.456 | 0.455
45 | 0.507 | 0.497 | 0.494 | 0.488 | 0.370 | 0.366 | 0.436 | 0.434
$\Delta\epsilon=0.25$
0 | 0.559 | 0.553 | 0.492 | 0.481 | 0.469 | 0.462 | 0.548 | 0.542
45 | 0.598 | 0.592 | 0.595 | 0.590 | 0.406 | 0.395 | 0.525 | 0.521
Table 6: Variation of critical buckling load with respect to the angle of
orientation of the elliptical cutout for different ply orientations for
$\gamma=10$ $\Delta\epsilon=0$
---
$\theta$ |  | Antisymmetric
---
Cross-Ply
| Symmetric
---
Cross-Ply
| Symmetric
---
Angle-Ply
| Antisymmetric
---
Angle-Ply
| LSM | MPM | LSM | MPM | LSM | MPM | LSM | MPM
0 | 0.488 | 0.485 | 0.417 | 0.409 | 0.443 | 0.443 | 0.469 | 0.468
45 | 0.509 | 0.509 | 0.502 | 0.495 | 0.83 | 0.378 | 0.448 | 0.449
$\Delta\epsilon=0.25$
0 | 0.590 | 0.581 | 0.509 | 0.501 | 0.485 | 0.479 | 0.576 | 0.572
45 | 0.619 | 0.612 | 0.618 | 0.610 | 0.410 | 0.405 | 0.560 | 0.554

###  3.6 Curvilinearly stiffened composite panels with a complicated cutout

The present method is applied to analyze a stiffened laminated composite plate
with a clover shaped cutout Figure 10. The cutout is a boolean of three
circles of radius $0.15$ resulting in a clover shape. The center of circles A,
B and C are at ($0.4,0.65$), ($0.5,0.7$) and ($0.5,0.6$) respectively. Figure
10 shows the mesh plot for the plate geometry with the red region being the
enriched elements.

To check the accuracy of the level set algorithm, results for clamped
curvilinearly stiffened composite plate with a clover shaped cutout are
compared with those obtained using ABAQUS. The curvilinear stiffener
configuration was adopted from [39] and scaled proportionally to match the
dimensions of the plate. From the ABAQUS library, the plate was meshed using
3040 S8R elements and curvilinear stiffener was meshed using 192 elements S8R
elements (48 divisions along the arc length and 4 divisions along the depth
direction). The IGA plate model was meshed using 4096 elements and the
curvilinear stiffener using 32 elements. The first five eigen buckling models
were extracted for each case. The buckling mode shapes obtained using the IGA
approach are similar to those obtained using ABAQUS as can be seen from Figure
11.

Figure 10: Isogeometric mesh plot of a plate with clover shaped cutout constructed with three cirles. The enriched elements are shown in red | Present | ABAQUS
---|---|---
|  |  |  |
Mode |  |  $\Delta T_{cr}$ |  |  $\Delta T_{cr}$
1 |  |  0.751 |  |  0.748
2 |  |  0.827 |  |  0.817
3 |  |  0.961 |  |  0.955
4 |  |  1.009 |  |  1
5 |  |  1.096 |  |  1.085
Figure 11: First five eigenmode shape plots for a plate with curvilinear
stiffener using a LSM and ABAQUS, a commercial available software.

##  4 Advantages and limitations

One of the advantages of the level set method is their ability to use the same
mesh controls on the plate for different types of cutouts. Figure 12 shows how
the plate is discretized using $16$ x $16$ elements for circular as well as
elliptical cutouts.

(a) circular

(b) angled elliptical

Figure 12: Identical IGA mesh for circular and elliptical cutouts

Discussions in Section 2.6 describe the creation of additional degrees of
freedom while implementing the level set method. Table 8 shows the number of
additional degrees of freedom created for each levels of refinement for a
plate with a cutout of radius $0.15$. It can be observed that the total
degrees of freedom (nDoF = Total number of Gauss points x 5 ) increase much
steeply than the additional degrees of freedom (Extra DoF) with increase in
refinement levels.

Table 8: Additional degrees of freedom with respect to an increase in refinement Refinement | nDoF | Extra DoF
---|---|---
3 | 660 | 160
4 | 1920 | 300
5 | 6320 | 540

Next, computational cost due to these additional degrees of freedom was
analyzed. Two cases were considered, a stiffened plate without any cutout
(Case 1) and a stiffened plate with a central circular cutout of radius
$0.15$. The material properties, ply orientation and stiffener profile were
the same as mentioned in Section 3.4. The computational time (in seconds) to
perform the static and the eigenvalue analyses are compared for each case
(Table 9, Figures 14 and 13). It is observed that the additional degrees of
freedom have a much higher influence on the results from a static analysis as
compared to the buckling analysis. However, the difference is very little for
lower levels of refinement.

Table 9: Computational time (in seconds) with respect to refinement for static and eigenvalue analyses Refinement | 4 | 5 | 6
---|---|---|---
| Static | Eigenvalue | Static | Eigenvalue | Static | Eigenvalue
Case 1 | 0.0342 | 0.2045 | 0.2008 | 0.7087 | 1.2722 | 4.3813
Case 2 | 0.0504 | 0.2244 | 0.2014 | 0.8724 | 4.5132 | 4.5352
$3$$4$$5$$6$$0$$2$$4$RefinementTime in secondsSimple platePlate with a cutout
Figure 13: Variation in computational time (in seconds) with respect to
refinement for eigenvalue analysis $3$$4$$5$$6$$0$$2$$4$RefinementTime in
secondsSimple platePlate with a cutout Figure 14: Variation in computational
time (in seconds) with respect to refinement for static analysis

##  5 Convergence studies

Convergence studies on a stiffened plate with a cutout are not as
straightforward as a simple plate problem since the size of the cutout, size
of the stiffener and the plate meshes all affect the results. We consider a
stiffened plate with a central circular cutout. The material properties, ply
orientation and stiffener profile were the same as mentioned in Section 3.4.
For convergence studies with respect to the increasing cutout size, an
unstiffened plate was considered. It can be observed from Figure 15 that the
rate of convergence becomes slower with an increase in the cutout radius.

$3$$4$$5$$6$$0.4$$0.6$$0.8$$1$RefinementCritical Buckling
LoadRadius0.10.150.20.25 Figure 15: Variation of critical buckling load with
respect to various plate mesh sizes for cutouts of increasing radii

Figures 16, 17 and 18 present the change of critical buckling loads with
respect to various stiffener and plate mesh sizes for cutouts of varying
radii. It is observed that for cutouts with relatively smaller radius ($0.15$
and $0.2$), for each refinement levels of the plate mesh, the results converge
for the $4^{th}$ order of refinement of the stiffener mesh. It is also
observed that the plate mesh density has a higher influence on the results
compared to the stiffener mesh density. However, as the radius of the cutout
is increased to $0.25$, the stiffener mesh density starts to impact the
convergence of results to a much higher degree (Figure 16).

$2$$4$$6$$8$$0.88$$0.9$$0.92$$0.94$Refinement (Stiffener)Critical Buckling
LoadRefinement (Plate)456 Figure 16: Variation of critical buckling load with
respect to various stiffener and plate mesh sizes for $radius=0.25$
$2$$4$$6$$8$$0.6$$0.7$$0.8$Refinement (Stiffener)Critical Buckling
LoadRefinement (Plate)3456 Figure 17: Variation of critical buckling load with
respect to various stiffener and plate mesh sizes for $radius=0.2$
$2$$4$$6$$8$$0.55$$0.6$Refinement (Stiffener)Critical Buckling LoadRefinement
(Plate)3456 Figure 18: Variation of critical buckling load with respect to
various stiffener and plate mesh sizes for $radius=0.15$

##  6 Conclusions

This paper discusses the level set method which requires very little
computational effort to obtain a mesh independent description of the
geometrical shape of a complicated cutout. The level sets discussed in this
paper are implemented in the form of a signed distance function which is both
computationally inexpensive and easy to implement. Such an implementation is
greatly useful since the tensor product nature of the NURBS basis functions
only allows for global refinement and therefore treating the trimmed objects
such as internal cutouts is rather difficult. However, implementing the level
set method obviates the use of trimmed NURBS surface to describe complex
geometrical features like the cutouts. The efficiency of this method for
thermal buckling analysis is described with the help of numerical examples for
stiffened laminated plates containing complicated cutouts of various shapes
and dimensions.

##  7 Acknowledgement

The authors would like to thank the Institute for Critical Technology and
Applied Science (ICTAS) at Virginia Tech for providing financial support which
made this research possible.

## References

## References

  * [1] T. J. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39-41) (2005) 4135–4195.
  * [2] T. T. Yu, S. Yin, T. Q. Bui, S. Hirose, A simple FSDT-based isogeometric analysis for geometrically nonlinear analysis of functionally graded plates, Finite Elements in Analysis and Design 96 (C) (2015) 1–10.
  * [3] S. De, K. Singh, B. Alanbay, R. K. Kapania, R. Aguero, Structural optimization of truck front-frame under multiple load cases, in: ASME International Mechanical Engineering Congress and Exposition, Vol. 52187, American Society of Mechanical Engineers, 2018, p. V013T05A039.
  * [4] S. De, K. Singh, J. Seo, R. K. Kapania, E. Ostergaard, N. Angelini, R. Aguero, Structural design and optimization of commercial vehicles chassis under multiple load cases and constraints, in: AIAA Scitech 2019 Forum, 2019, p. 0705\.
  * [5] S. De, K. Singh, J. Seo, R. K. Kapania, E. Ostergaard, N. Angelini, R. Aguero, Lightweight chassis design of hybrid trucks considering multiple road conditions and constraints, World Electric Vehicle Journal 12 (1) (2021) 3.
  * [6] S. De, K. Singh, J. Seo, R. Kapania, R. Aguero, E. Ostergaard, N. Angelini, Unconventional truck chassis design with multi-functional cross members, Tech. rep., SAE Technical Paper (2019).
  * [7] M. Jrad, S. De, R. K. Kapania, Global-local aeroelastic optimization of internal structure of transport aircraft wing, in: 18th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, 2017, p. 4321.
  * [8] S. De, R. K. Kapania, Algorithms for 2d mesh decomposition in distributed design optimization (2021).  arXiv:2002.00525.
  * [9] Y. Bazilevs, V. M. Calo, T. J. R. Hughes, Y. Zhang, Isogeometric fluid-structure interaction: theory, algorithms, and computations, Computational Mechanics 43 (1) (2008) 3–37.
  * [10] J. A. Evans, Y. Bazilevs, I. Babuška, T. J. Hughes, $n$-Widths, $sup$-infs, and optimality ratios for the $k$-version of the isogeometric finite element method, Computer Methods in Applied Mechanics and Engineering 198 (21-26) (2009) 1726–1741.
  * [11] T. J. R. Hughes, A. Reali, G. Sangalli, Efficient quadrature for NURBS-based isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 199 (5-8) (2010) 301–313.
  * [12] F. Auricchio, F. Calabrò, T. J. Hughes, A. Reali, G. Sangalli, A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 249-252 (2012) 15–27.
  * [13] Y. Bazilevs, C. Michler, V. M. Calo, T. J. Hughes, Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly enforced boundary conditions on unstretched meshes, Computer Methods in Applied Mechanics and Engineering 199 (13-16) (2010) 780–790.
  * [14] N. Valizadeh, T. Bui, V. Vu, H. Thai, M. Nguyen, Isogeometric Simulation for Buckling, Free and Forced Vibration of Orthotropic Plates, International Journal of Applied Mechanics 05 (02) (2013) 1350017.
  * [15] S. Yin, J. S. Hale, T. Yu, T. Q. Bui, S. P. Bordas, Isogeometric locking-free plate element: A simple first order shear deformation theory for functionally graded plates, Composite Structures 118 (1) (2014) 121–138.
  * [16] H. Kapoor, R. K. Kapania, S. R. Soni, Interlaminar stress calculation in composite and sandwich plates in NURBS Isogeometric finite element analysis, Composite Structures 106 (2013) 537–548.
  * [17] H. Kapoor, R. K. Kapania, Geometrically nonlinear NURBS isogeometric finite element analysis of laminated composite plates, Composite Structures 94 (12) (2012) 3434–3447.
  * [18] C. V. Verhoosel, M. A. Scott, T. J. R. Hughes, R. de Borst, An isogeometric analysis approach to gradient damage models, International Journal for Numerical Methods in Engineering 86 (1) (2011) 115–134.  arXiv:1010.1724.
  * [19] J. Lu, Isogeometric contact analysis: Geometric basis and formulation for frictionless contact, Computer Methods in Applied Mechanics and Engineering 200 (5-8) (2011) 726–741.
  * [20] H. Nguyen-Xuan, L. V. Tran, C. H. Thai, S. Kulasegaram, S. P. Bordas, Isogeometric analysis of functionally graded plates using a refined plate theory, Composites Part B: Engineering 64 (2014) 222–234.
  * [21] W. A. Wall, M. A. Frenzel, C. Cyron, Isogeometric structural shape optimization, Computer Methods in Applied Mechanics and Engineering 197 (33-40) (2008) 2976–2988.
  * [22] J. Kiendl, K.-U. Bletzinger, J. Linhard, R. Wüchner, Isogeometric shell analysis with Kirchhoff–Love elements, Computer Methods in Applied Mechanics and Engineering 198 (49-52) (2009) 3902–3914.
  * [23] J. Kiendl, Y. Bazilevs, M. C. Hsu, R. Wüchner, K. U. Bletzinger, The bending strip method for isogeometric analysis of Kirchhoff-Love shell structures comprised of multiple patches, Computer Methods in Applied Mechanics and Engineering 199 (37-40) (2010) 2403–2416.
  * [24] S. Shojaee, E. Izadpanah, N. Valizadeh, J. Kiendl, Free vibration analysis of thin plates by using a NURBS-based isogeometric approach, Finite Elements in Analysis and Design 61 (2012) 23–34.
  * [25] L. V. Tran, C. H. Thai, H. T. Le, B. S. Gan, J. Lee, H. Nguyen-Xuan, Isogeometric analysis of laminated composite plates based on a four-variable refined plate theory, Engineering Analysis with Boundary Elements 47 (2014) 68–81.
  * [26] L. A. Piegl, W. Tiller, Geometry-based triangulation of trimmed NURBS surfaces, Computer-Aided Design 30 (1) (1998) 11–18.
  * [27] S. Shojaee, N. Valizadeh, E. Izadpanah, T. Bui, T. V. Vu, Free vibration and buckling analysis of laminated composite plates using the NURBS-based isogeometric finite element method, Composite Structures 94 (5) (2012) 1677–1693.
  * [28] R. Schmidt, R. Wüchner, K. U. Bletzinger, Isogeometric analysis of trimmed NURBS geometries, Computer Methods in Applied Mechanics and Engineering 241-244 (2012) 93–111.
  * [29] S. S. Ghorashi, N. Valizadeh, S. Mohammadi, Extended isogeometric analysis for simulation of stationary and propagating cracks, International Journal for Numerical Methods in Engineering 89 (9) (2012) 1069–1101.  arXiv:1010.1724.
  * [30] E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, A. Düster, Geometric modeling, isogeometric analysis and the finite cell method, Computer Methods in Applied Mechanics and Engineering 249-252 (2012) 104–115.
  * [31] D. Schillinger, L. Dedè, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, T. J. Hughes, An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Computer Methods in Applied Mechanics and Engineering 249-252 (2012) 116–150.
  * [32] B. Devarajan, D. Locatelli, R. Kapania, R. Meritt, Thermo-mechanical analysis and design of threaded fasteners, in: 2016 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2016, p. 0579.
  * [33] Q. Li, B. Devarajan, X. Zhang, R. Burgos, D. Boroyevich, P. Raj, Conceptual design and weight optimization of aircraft power systems with high-peak pulsed power loads, in: SAE Technical Paper, SAE International, 2016, p. 10.  doi:10.4271/2016-01-1986.
  * [34] A. K. Noor, W. S. Burton, C. W. Bert, Computational Models for Sandwich Panels and Shells, Applied Mechanics Reviews 49 (3) (1996) 155.
  * [35] W. S. Burton, A. K. Noor, Three Dimensional Solutions for Thermomechanical Stresses in Sandwich Panels and Shells, Journal of Engineering Mechanics 120 (10) (1994) 2044–2071.
  * [36] A. K. Noor, W. S. Burton, Three Dimensional Solutions for Thermal Buckling of Multilayered Anisotropic Plates, Journal of Engineering Mechanics 118 (4) (1992) 683–701.
  * [37] M. Mukhopadhyay, A. Mukherjee, Finite element buckling analysis of stiffened plates, Computers & structures 34 (6) (1990) 795–803.
  * [38] B. G. Prusty, S. Satsangi, Finite element transient dynamic analysis of laminated stiffened shells, Journal of Sound and Vibration 248 (2) (2001) 215–233.
  * [39] A. Y. Tamijani, R. K. Kapania, Buckling and Static Analysis of Curvilinearly Stiffened Plates Using Mesh-Free Method, AIAA Journal 48 (12) (2010) 2739–2751.
  * [40] W. Zhao, R. K. Kapania, Buckling analysis of unitized curvilinearly stiffened composite panels, Composite Structures 135 (2016) 365–382.
  * [41] J. H. Robinson, S. Doyle, G. Ogawa, M. Baker, S. De, M. Jrad, R. K. Kapania, Aeroelastic optimization of wing structure using curvilinear spars and ribs (sparibs), in: 17th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, 2016, p. 3994.
  * [42] S. De, M. Jrad, D. Locatelli, R. K. Kapania, M. Baker, Sparibs geometry parameterization for wings with multiple sections using single design space, in: 58th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2017, p. 0570.
  * [43] S. De, Structural modeling and optimization of aircraft wings having curvilinear spars and ribs (sparibs), Ph.D. thesis, Virginia Tech (2017).
  * [44] P. Shi, R. K. Kapania, C. Y. Dong, Vibration and Buckling Analysis of Curvilinearly Stiffened Plates Using Finite Element Method, AIAA Journal 53 (5) (2015) 1319–1335.
  * [45] P. Hao, B. Wang, K. Tian, G. Li, X. Zhang, Optimization of curvilinearly stiffened panels with single cutout concerning the collapse load, International Journal of Structural Stability and Dynamics 16 (07) (2016) 1550036\.
  * [46] B. Devarajan, R. K. Kapania, Thermal buckling of curvilinearly stiffened laminated composite plates with cutouts using isogeometric analysis, Composite Structures (2020) 111881.
  * [47] B. Devarajan, Thermomechanical and Vibration Analysis of Stiffened Unitized Structures and Threaded Fasteners, vtechworks.lib.vt.edu, 2019.
  * [48] L. L. Piegl, W. Tiller, The NURBS Book, Vol. 28, Springer Science & Business Media, 1996.  arXiv:arXiv:1011.1669v3.
  * [49] L. V. Tran, M. A. Wahab, S. E. Kim, An isogeometric finite element approach for thermal bending and buckling analyses of laminated composite plates, Composite Structures 179 (2017) 35–49.
  * [50] C. Sarath Babu, T. Kant, Refined higher order finite element models for thermal buckling of laminated composite and sandwich plates, Journal of Thermal Stresses 23 (2) (2000) 111–130.
  * [51] T. Kant, C. S. Babu, Thermal buckling analysis of skew fibre-reinforced composite and sandwich plates using shear deformable finite element models, Composite Structures 49 (1) (2000) 77–85.
  * [52] A. Avci, S. Kaya, B. Daghan, Thermal buckling of rectangular laminated plates with a hole, Journal of Reinforced Plastics and Composites 24 (3) (2005) 259–272.