A cutout isogeometric analysis for thin laminated composite plates using level sets

Composite Structures 127 (2015) 152–164

Contents lists available at ScienceDirect

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

A cutout isogeometric analysis for thin laminated composite plates using level sets Shuohui Yin a, Tiantang Yu a,⇑, Tinh Quoc Bui b,*, Shifeng Xia a,c, Sohichi Hirose b a

Department of Engineering Mechanics, Hohai University, Nanjing 210098, PR China Department of Mechanical and Environmental Informatics, Graduate School of Information Science and Engineering, Tokyo Institute of Technology, 2-12-1-W8-22, Ookayama, Meguro-ku, Tokyo 152-8552, Japan c College of Energy and Electrical Engineering, Hohai University, Nanjing 210098, PR China b

a r t i c l e

i n f o

Article history: Available online 14 March 2015 Keywords: Isogeometric analysis NURBS Laminates Plates Free vibration Buckling

a b s t r a c t Numerical modeling with treatment of trimmed objects such as internal cutouts in terms of NURBSbased isogeometric analysis presents several challenges, primarily due to the tensor product of the NURBS basis functions. In this paper we develop a new simple and effective isogeometric analysis for modeling buckling and free vibration problems of thin laminated composite plates with cutouts. We adopt the classical plate theory for the present formulation. The new approach can nicely 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. Numerical examples with complicated shapes are considered and analyzed to show the influences of cutout geometry, fiber orientation, boundary conditions, etc. on natural frequency and buckling behaviors of laminated plates. The results are compared with reference solutions showing a high accuracy of the proposed method. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Owing to the advantages like light weight, longer life, and fatigue endurance etc., thin laminated composites have been widely used as components in engineering structures. Due to a large requirement of laminates for a variety of engineering applications, thin plate structures with arbitrary cutouts are inevitable. The presence of cutouts can significantly affect the dynamic and buckling behaviors of structures in particular, and their performance in general. Studies on the analysis of eigenvalues and stability of thin laminates with cutout are of great importance for many practical applications including airplane wing, fuselage and ribs. Because of the complexity of the laminated composites with cutouts, numerical methods are extensively used in this subject. Studies on the vibration and buckling problems of thin laminated composite plates with cutouts have performed using different numerical approaches including finite element method (FEM) [1–4], Rayleigh–Ritz method [5], meshfree method [6], finite strip method [7,8] and extended finite element method (XFEM) [9]. In ⇑ Corresponding authors. Tel.: +86 2552430342 (T.T. Yu), +81 357343587 (T.Q. Bui). E-mail addresses: [email protected] (T. Yu), [email protected] (T.Q. Bui). http://dx.doi.org/10.1016/j.compstruct.2015.03.016 0263-8223/Ó 2015 Elsevier Ltd. All rights reserved.

recent years, the isogeometric analysis (IGA) [10,11] has introduced and 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. Some insights into mathematical properties [12,13], integration method [14,15] and splines techniques [16] have been studied. The inherent characteristics of IGA make it superior to the classical FEM in many aspects [10]. The IGA has successfully applied to many engineering problems including plates and shells [17–20], fluid mechanics [21], fluid– structure interaction [22], damage and fracture mechanics [23], contact mechanics [24], unsaturated flow problem in porous media [25], and structural shape optimization [26]. Interesting features of exact geometry representation and highorder smoothness of NURBS basis functions have attractive for the analysis of plates and shells problems. The C1 continuity requirement of classical plate theory is easily handled without additional efforts due to the high-order continuity of NURBS basis functions. Based on classical plate theory, a direct construction of rotationfree isogeometric shell was initially introduced in [27], and fully developed for multiple NURBS patches using the bending strip method by Kiendl [28]. Later, the rotation-free isogeometric shell element was extended to large deformation [29], free vibration and buckling analysis of laminated plates [30] and functional

153

S. Yin et al. / Composite Structures 127 (2015) 152–164

graded plates [31], cloth simulation [32]. Investigated in the previous works, it exhibits that the rotation-free isogeometric plate/ shell elements can attain very good accuracy and are efficient for plate and shell structures. However, most of the applications are limited to single-patch structures and simply geometry. For complex geometries like structure with cutouts, the trimmed NURBS surface is useful. Shojaee et al. [30,33] artificially divided the square plate with a hole of complicated shape into several NURBS patches and then applied the bending strip method to maintain C1 continuity between patches. This approach however is found to be ineffective and not straightforward to unify design and analysis. To deal with the trimmed NURBS surfaces, Schmidt et al. [34] alternatively developed a local reconstruction technique to rebuild trimmed elements with separate patches in which the geometric errors and loss of the higher-order continuity along the interior and limited to single trimming curves are introduced. Nevertheless, numerical modeling with treatment of trimmed objects such as internal cutouts in terms of NURBS-based IGA has some drawbacks primarily because of the tensor product induced by the NURBS basis. The existing geometric models with trimmed NURBS surface thus remain a challenging problem in computational mechanics community. In order to conveniently deal with complex geometry problems, the IGA on one hand has combined with the enrichment method to form the so-called XIGA [35], and on the other hand that combines with the finite cell method (FCM) resulting the finite cell extension to isogeometric analysis [36,37], in which the level sets are used to define the geometry of the computational domain. In this work, the level sets are used to describe the internal cutouts and the numerical integration is used only inside the physical domain. One must be noted that most of the preceding studies using IGA with level sets were developed for elastoelastic problems, but have not accounted yet for dynamic and buckling problems of thin laminated composite structures with complicated shapes. The objective of this paper is to fill such a gap which essentially devotes to modeling dynamic and buckling of complicated thin laminated structures by a new cutout IGA approach. This new simple and effective IGA is found accurately and effectively, which will be demonstrated later in the numerical results, in modeling complicated thin laminated composites with cutouts based on the classical plate theory. The rest of this paper is structured as follows: fundamental of NURBS basis functions and their derivatives are presented in Section 2. In Section 3, the IGA discrete equations based on classical plate theory and level sets for thin laminated composites with cutouts are derived. Numerical examples with complicated shapes for the free vibration problems of thin laminated structures with cutouts are presented in the next section to show the accuracy of the proposed approach. Similarly, buckling results for thin laminates with cutouts are presented in Section 5. We end with some conclusions drawn from the study in Section 6.

With a given knot vector kðnÞ, the B-spline basis function, written as N i;p ðnÞ, is defined recursively as follows [38]:

Ni;0 ðnÞ ¼



1 if ni 6 n < niþ1 0

otherwise

for p ¼ 0;

and

Ni;p ðnÞ ¼

niþpþ1  n n  ni Ni;p1 ðnÞ þ Niþ1;p1 ðnÞ for p P 1: niþp  ni niþpþ1  niþ1 ð2Þ

A quadratic B-spline basis functions example is presented in Fig. 1 using the following knot vector

kðnÞ ¼ fn1 ; n2 ; . . . n9 ; n10 gT ¼ f0; 0; 0; 0:2; 0:4; 0:6; 0:8; 1; 1; 1gT It can be observed that a B-Spline basis function is C

1

ð3Þ

continu-

ous inside a knot span, i.e., between two distinct knots, and C p1 continuous at a single knot. This character satisfies the requirement of the classical plate theory. The NURBS basis function Ri;p ðnÞ in the framework of partition of unity is constructed by a weighted average of the B-spline basis functions [38] as follows:

Ni;p ðnÞwi Ri;p ðnÞ ¼ Pn j¼1 N j;p ðnÞwj

ð4Þ

where wi is the ith weight; the NURBS basis function is degenerate into B-spline basis function for wi ¼ 1. Similarly, the bivariate NURBS basis function for a NURBS surface is given by

Ni;p ðnÞNj;q ðgÞwi;j Ni;p ðnÞN j;q ðgÞwi;j Rp;q ¼ i;j ðn; gÞ ¼ Pn Pm Wðn; gÞ i¼1 j¼1 N i;p ðnÞN j;q ðgÞwi;j

ð5Þ

where wi;j represents the 2D weight; Nj;q ðgÞ is the B-spline basis of order p defined on the knot vector kðgÞ, followed the recursive forP P mula shown in Eqs. (1) and (2); Wðxi; gÞ ¼ ni¼1 m j¼1 N i;p ðnÞN j;q ðgÞwi;j is the weighting function for a NURBS surface. It should be stated here that the NURBS basis functions own the same properties as B-splines. By using the NURBS basis functions, a NURBS surface of order p in the n direction and order q in the g direction can be constructed as follows:

Sðn; gÞ ¼

n X m X Rp;q i;j ðn; gÞP i;j ;

ð6Þ

i¼1 j¼1

where P i;j represents the coordinates of control points in two dimensions.

2. NURBS basis functions and their derivatives The NURBS basis functions and their derivatives are briefly introduced in this section. For details, curious readers may refer to [10,38]. The NURBS is known as a generalization of B-splines. A B-spline is defined on a one dimensional parametric space n 2 ½0; 1, by a set of non-decreasing numbers called knot vector kðnÞ ¼ fn1 ¼ 0; . . . ; ni ; . . . ; nnþpþ1 ¼ 1gT ðni 6 niþ1 Þ, and a set of control points P i ; i ¼ 1; . . . ; n, where n and p is the number of splines basis functions and the order of splines basis functions, respectively. The non-zeros knot span ½ni ; niþ1 Þ is defined as element in IGA. A knot vector kðnÞ is called an open knot vector when the two ends of the knot are repeated p þ 1 times.

ð1Þ

Fig. 1. Quadratic B-spline basis functions.

154

S. Yin et al. / Composite Structures 127 (2015) 152–164

The first derivatives of the NURBS basis function Rp;q i;j ðxi; gÞ with respect to each parametric variable are derived by applying the quotient rule to Eq. (5) as @N i;p ðnÞ Nj;q ð @n

gÞ gÞWðxi; gÞ  @Wðxi; Ni;p ðnÞN j;q ðgÞ @n 2 @n ðWðxi; gÞÞ @N j;q ðgÞ gÞ p;q Ni;p ðnÞWðxi; gÞ  @Wðxi; Ni;p ðnÞNj;q ðgÞ @Ri;j ðxi; gÞ @g @g ¼ wi;j 2 @g ðWðxi; gÞÞ

@Rp;q i;j ðxi; gÞ

¼ wi;j

ð7aÞ ð7bÞ

The first derivatives of the weighting function Wðxi; gÞ with respect to each parametric variable are given by n X m @Wðxi; gÞ X @Ni;p ðnÞ ¼ N j;q ðgÞwi;j @n @n i¼1 j¼1

ð8aÞ

n X m @Wðxi; gÞ X @Nj;q ðgÞ ¼ Ni;p ðnÞwi;j @g @n i¼1 j¼1

ð8bÞ

Due to the assumption of classical plate theory for a thin laminated composite plate, D can be written as

2

D11

6 D ¼ 4 D12 D16

D12

D16

3

D22

7 D26 5

D26

D66

ð14Þ

with

DIJ ¼

N   1X ðQ IJ Þk z3k  z3k1 ; 3 k¼1

I; J ¼ 1; 2; 6

ð15Þ

where N is the number of layers of the composite laminated plate  IJ are obtained by and Q 2

4

Q 11 ¼ Q 11 cos4 h þ 2ðQ 12 þ 2Q 66 Þ sin h cos2 h þ Q 22 sin h 2

ð16aÞ

4

Q 12 ¼ ðQ 11 þ Q 22  4Q 66 Þ sin h cos2 h þ Q 12 ðsin h þ cos4 hÞ

Higher-order derivatives of the NURBS basis function and the weighting function can be obtained in a similar process.

ð16bÞ 3

Q 16 ¼ ðQ 11  Q 12  2Q 66 Þ sin h cos3 h þ ðQ 12  Q 22 þ 2Q 66 Þ sin h cos h 4

2

Q 22 ¼ Q 11 sin h þ 2ðQ 12 þ 2Q 66 Þ sin h cos2 h þ Q 22 cos4 h

3. Isogeometric analysis for thin laminated composite plates with cutouts

Q 26 ¼ ðQ 11  Q 12  2Q 66 Þ sin h cos h þ ðQ 12  Q 22 þ 2Q 66 Þ sin h cos3 h

3.1. Governing equations and discretization

Q 66 ¼ ðQ 11 þ Q 22  2Q 12  2Q 66 Þ sin h cos2 h þ Q 66 ðsin h þ cos4 hÞ

For convenience, the control points are renumbered with a unified subscript I based on the subscripts {i; j} in Eq. (6), and the NURBS surface are represented as

with

xðnÞ ¼

NP X RI ðnÞP I

3

2

ð9Þ

NP X wðxðnÞÞ ¼ RI ðnÞwI

ð10Þ

Q 12 ¼

Q 22

Q 66 ¼ G12 ;

Z X

deTp rp dX ¼

8 9 @w > < z @x > = u¼ v ¼ z @w @y > > : > ; > : ; w w

ð11Þ

and the pseudo-strains and pseudo-stresses are

ep

8 @2 w 9  2 > > > < @x2 > = ¼  @@yw2 ; > > > > : @2 w ; 2 @x@y

8 9 > < Mx > = rp ¼ M y > > : ; M xy

ð12Þ

The relationship between the pseudo-strains and pseudostresses are derived from Hooke’s law by

rp ¼ Dep

ð13Þ

ð16fÞ

v 21 E1 ¼ v 12 E2

ð17Þ

Z

€ dX duT qu

ð18Þ

X

In the case of in-plane buckling analyses and assuming prebuckling stresses r0 , the weak-form can be briefly expressed as

Z X

deTp rp dX þ h 3

8 9 > =

ð16eÞ

in which, h is the fiber orientation of each layer; E1 and E2 are the Young modulus in the x and y directions, respectively; v 12 and v 21 are the Poisson’s ratios; and G12 is the shear modulus. The weak-form of dynamic equations of the laminated composite plate without external forces can be derived from the principle of virtual work as

I¼1

where wI is the deflection of control point I. The displacements of an arbitrary point in plate are then expressed as

ð16dÞ

v 12 E2 ; 1  v 12 v 21

Q 11 ¼

I¼1

where x ¼ ðx; yÞ is the physical coordinates vector; n ¼ ðxi; gÞ is the parametric coordinates vector; NP represents the total number of control points; RI ðnÞ is the NURBS basis function defined in Eq. (5) and the superscripts of basis function are ignored. As higher order basis functions can be easily obtained in IGA, the relationship between deflection and rotation can be adopted directly. By employing the classical plate theory, the deflection of the plate wðxÞ is used as the only independent variable and the other two displacement components can be obtained through the derivatives of wðxÞ. The approximation of the deflection field is given by

E1 ; 1  v 12 v 21 E2 ¼ ; 1  v 12 v 21

4

ð16cÞ

þ

h 12

Z

rT d

X

Z

3

rT dwr0 rwdX þ

X

h 12

Z X

rT d

@w @w r0 r dX @x @x

@w @w r0 r dX ¼ 0 @y @y

ð19Þ "

#

r0x s0xy in which r ¼ ½ @=@x @=@y  and r0 ¼ are the gradient s0xy r0y T

operator and in-plane pre-buckling stresses, respectively, and h is plate thickness. Substituting the displacements u, pseudo-strains ep and pseudo-stresses rp from Eqs. (12), (13) and (15) into Eqs. (25) and (26), the eigenvalue equations of the free vibration and buckling analysis are written in the following forms

ðK  x2 MÞU ¼ 0

ð20Þ

ðK  Ncr K G ÞU ¼ 0

ð21Þ

where x is the natural frequency, Ncr is the critical buckling load, U is the eigenvector fw1 ; w2 ; . . . ; wNP g consisted of deflection at all control points, K; M and K G are the global stiffness, mass and geometrical stiffness matrix given by

155

S. Yin et al. / Composite Structures 127 (2015) 152–164

Fig. 2. The process of geometric modeling for a geometric model with a complex cutout.

K¼

Z

BT DBdX

ð22Þ

~ T qBd ~ X B

ð23Þ

X

M¼

Z

X G

K ¼h

Z

3

X

GT1

h r0 G1 dX þ 12

Z

in which, at control point I

X

3

GT2

h r0 G2 dX þ 12

Z X

BI ¼ f RI;xx

RI;yy

2RI;xy gT

ð25Þ

~ I ¼ f zRI;x B

zRI;y

RI gT

ð26Þ

G1I ¼ f RI;x GT3

r0 G3 dX

ð24Þ

G2I ¼ f RI;xx G3I ¼ f RI;xy

T

RI;y g

RI;xy g

ð27Þ T

ð28Þ

T

ð29Þ

RI;yy g

156

S. Yin et al. / Composite Structures 127 (2015) 152–164

Fig. 3. A square plate with a single circular hole.

3.2. Level sets for complex cutout and its implementation A geometric model with complex cutout is tackled via Boolean operations, i.e. union, intersection and difference, it hence results in a trimmed NURBS object, see Fig. 2 for a schematic representation of the object. As mentioned in [34] trimmed NURBS object in the CAD system are only visually trimmed, the parameterization remains unchanged and the cropped part is not visualized. That implies that the geometric information of the structural object has changed, but the NURBS topology has not altered and the entire construction process is being stored. The method described in [34] for trimmed NURBS surface is in general not effective and limited to single trimming curves. Fig. 2 shows a construction process of geometric modeling for a geometric model with a complex cutout. Based on the initial geometry X, the cutout is composed by several interacted trimming curves, for instance, C1 ; C2 and C3 . We can see in the picture that each Boolean operation corresponds to a trimming curve (the rectangle can be viewed as a curve consisting of four connected straight lines). On the other hand, each trimming curve can be

(a) 15 ×15 control points

described by a level set ui such as ðC1 ! u1 ; C2 ! u2 ; C3 ! u3 Þ. As the complex cutout can be implicitly described by the process of Boolean operation, we can implicitly transfer the process of Boolean operation on the level sets that corresponds to each trimming curve. That is to say that each Boolean operation on the NURBS geometry corresponding to a level set, and the entire Boolean operation on the NURBS geometry can be transferred to the operation on the level sets. Thus, the final geometry for a geometric model with a complex cutout can be described as Xfinal ¼ X=ðu1 [ u2 [ u3 Þ and the final level set to describe the cutout and differentiate between inside and outside the physical domain is u ¼ u1 [ u2 [ u3 . For the geometric model with complex cutout shown in Fig. 2, the computational mesh for IGA is generated with a square domain, not considering the cutout yet. In the analysis, the location of an integration point xp is judged according to the three level sets, i.e., u1 ðxp Þ; u2 ðxp Þ, and u3 ðxp Þ, integration point xp locates inside the physical domain when u1 ðxp Þ < 0 or u2 ðxp Þ < 0 or u3 ðxp Þ < 0, while integration point xp locates outside the physical domain when u1 ðxp Þ > 0 and u2 ðxp Þ > 0 and u3 ðxp Þ > 0. Based on the basic concepts fictitious domain approaches [36,37], the level sets are used to differentiate between inside and outside the physical domain and the numerical integration is used only inside the physical domain. In order to improve the accuracy of integration in the element split by cutout, some integration approaches are developed such as triangular sub-domain technique used in XIGA [35,39] and XFEM [40,41], composed adaptive integration used in FCM [36,37] and the new numerical integration proposed in [42]. In this paper, the triangular sub-domain technique however is adopted. Consequently, the weak-form of the free vibration and buckling analysis Eqs. (25) and (26) can be expressed as

Z X

Z

X

deTp arp dX ¼

deTp arp dX þ h 3

þ

Z

h 12

Z

€ dX duT aqu

ð30Þ

X

Z

rT dwar0 rwdX

X 3

rT d

X

@w @w h ar0 r dX þ @x @x 12

Z X

rT d

@w @w ar0 r dX ¼ 0 @y @y

ð31Þ

where aðxÞ is set to a very small value (1012 in this paper) when the integral point x is inside of the cutout, aðxÞ ¼ 1 when the integral point x is inside of physical domain.

(b) 27 × 27 control points

Fig. 4. NURBS meshes for the square plate, green circles denote control points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

157

S. Yin et al. / Composite Structures 127 (2015) 152–164 Table 1 Normalized frequencies of clamped square plate with a circular hole. Mode

This work Number of control points

1 2 3 4 5 6 7 8

99

15  15

21  21

27  27

33  33

6.237 8.796 8.796 10.552 11.618 12.158 13.179 13.179

6.186 8.748 8.748 10.521 11.567 12.066 13.095 13.095

6.181 8.675 8.675 10.517 11.560 12.045 13.036 13.036

6.184 8.657 8.657 10.513 11.557 12.055 13.021 13.021

6.183 8.657 8.657 10.513 11.559 12.052 13.022 13.022

Table 2 Natural frequencies of a completely free square plate with four holes.

Ref. [43]

Ref. [6]

6.24 8.457 8.462 10.233 11.719 12.299 13.037 13.041

6.149 8.577 8.634 10.422 11.414 11.838 12.829 12.842

Based on the above concept, the eigenvalue equations can be obtained as Eqs. (20) and (21), and the global stiffness K matrix, mass matrix M and geometrical stiffness matrix K G are reformulated as

Z K¼ BT aDBdX ZX ~ T aqBd ~ X M¼ B

ð32Þ ð33Þ

X

Z 3Z 3Z h h GT2 ar0 G2 dXþ GT ar0 G3 dX K G ¼h GT1 ar0 G1 dXþ 12 X 12 X 3 X

ð34Þ

Mode

ABAQUS

This work

4 5 6 7 8 9 10

6.036 8.586 10.070 15.651 15.651 26.693 26.693

6.027 8.573 10.026 15.620 15.622 26.608 26.609

Table 3 Natural frequencies of a simply supported square plate with four holes. Mode

ABAQUS

This work

1 2 3 4 5 6 7 8 9 10

9.016 22.704 22.705 38.493 43.778 44.275 57.234 57.238 75.691 75.694

8.997 22.652 22.653 38.468 43.552 44.150 56.930 56.931 75.490 75.493

Table 4 Natural frequencies of a clamped supported square plate with four holes.

4. Free vibration analysis The verification and the accuracy of the present approach for free vibration analysis of different complicated geometrical shapes with various boundary conditions are investigated. Throughout the study, cubic NURBS basis function is adopted. For p-order basis functions, p + 1 Gauss points per element are needed in order to integrate both mass and stiffness matrices accurately [14]. A 4  4 Gaussian quadrature scheme is hence used for ordinary elements; while a triangular sub-domain technique is adopted in the elements cut by internal cutout to improve the accuracy, and a 5order Gaussian quadrature scheme is used in each triangular subdomain. Different boundaries of the plate are named as simply supported (S), clamped (C), and free (F). For the clamped boundary conditions, the rotations are obtained from the derivatives of transverse deflection. Consequently, the constraint on the rotations is hence imposed by fixing the transverse defection with two rows of control points as described in [30].

(a) Geometry (unit:m)

Mode

ABAQUS

This work

1 2 3 4 5 6 7 8 9 10

17.504 36.032 36.034 58.239 58.390 63.796 77.246 77.251 93.175 93.179

17.492 35.987 35.988 57.938 58.547 63.785 77.012 77.026 92.869 92.871

The typical boundary conditions can be rewritten as follows: Simply supported edges:

w¼0

ð35Þ

(b) Mesh

Fig. 5. The geometry and mesh of a square plate with four holes.

158

S. Yin et al. / Composite Structures 127 (2015) 152–164

1th Mode

2th Mode

3th Mode

4th Mode

5th Mode

6th Mode

Fig. 6. First six mode shapes of a CCCC square plate with four holes.

Table 5 Normalized frequencies of a completely free square plate with a complicated cutout. Mode

IGA with 8-patches [33]

MKI [44]

NS-RPIM [45]

This work

4 5 6 7 8 9 10 11 12 13

3.482 3.968 5.216 5.805 6.015 7.100 7.485 8.114 8.853 9.001

3.393 4.277 5.936 7.639 8.005 9.225 9.893 – – –

3.234 3.846 4.684 5.497 5.547 6.929 7.350 8.103 8.392 8.667

3.223 3.823 4.657 5.472 5.527 6.924 7.336 8.076 8.357 8.664

the radius is R = 1 m, Young’s modulus E ¼ 2:0  1011 N/m2, Poisson’s ratio m ¼ 0:3, density q ¼ 8000 kg/m2. The fully clamped boundary condition is considered and five different meshes with 9  9; 15  15; 21  21; 27  27 and 33  33 control points are studied. Fig. 4 shows two different meshes of 15  15 and 27  27 control points. For convenience in representing the results, the frequency ~ ¼ ½qhx2 a4 =Dð1  v 2 Þ1=4 with parameter is normalized by x Fig. 7. Geometric parameters of a square plate with a complicated cutout.

Clamped edges:

w ¼ @w=@x ¼ @w=@y ¼ 0

ð36Þ

4.1. Isotropic plates with single and multiple circular holes We start by considering an isotropic square plate with a single circular hole as shown in Fig. 3. The geometry and material parameters are taken as follows: a = 10 m, thickness h = 0.1 m,

3

D ¼ Eh =12ð1  v 2 Þ. The first eight frequencies computed by the present method are tabulated in Table 1 and are compared with the reference solutions in [6,43]. As expected, a good convergence and agreement among them is found. Next example considers a square thin plate with four holes as depicted in Fig. 5 using 21  21 control points and 19  19 elements. The same material parameters are taken but the geometry parameters are as follows: a = 5 m, thickness h = 0.05 m, R = 0.5 m. The completely free, the simply and clamped supported boundary conditions are considered. For comparison purpose, a reference solution is obtained by using the finite element software ABAQUS

159

S. Yin et al. / Composite Structures 127 (2015) 152–164

with a very fine mesh (e.g., 2401 nodes, 2248 elements). Tables 2–4 show a comparison of natural frequencies calculated by the present method and by the FEM (ABAQUS). It is apparent that a good agreement between two methods is found. The first six eigenmodes of the clamped boundary conditions are also given in Fig. 6.

Table 6 Normalized frequencies of a simply supported square plate with a complicated cutout. Mode

IGA with 8-patches [33]

MKI [44]

NS-RPIM [43]

EFG [46]

This work

1 2 3 4 5 6 7 8 9 10

5.193 6.579 6.597 7.819 8.812 9.420 10.742 10.776 11.919 13.200

5.390 7.502 8.347 10.636 11.048 12.894 13.710 14.062 16.649 17.364

4.919 6.398 6.775 8.613 9.016 10.738 10.930 11.601 12.903 13.283

5.453 8.069 9.554 10.099 11.328 12.765 13.685 14.305 15.721 17.079

4.912 6.396 6.770 8.561 8.992 10.670 10.888 11.590 12.806 13.180

4.2. Isotropic plate with a complicated cutout The applicability of the present method to complicated cutout plate is analyzed by considering a thin isotropic plate with a complicated cutout. The geometry of the plate and its dimensions are shown in Fig. 7. Plate thickness is h = 0.05 m and the material parameters are the same as those considered in the previous subsection. The general process of geometric modeling of this example is presented in Fig. 2. The same example is analyzed by Shojaee et al. [33] using IGA and Kirchhoff–Love theory with eight NURBS patches and 512 control points, and the bending strip method is

Table 7 Normalized frequencies of a clamped square plate with a complicated cutout. Mode

IGA with 8-patches [33]

NS-RPIM [45]

EFG [46]

This work

1 2 3 4 5 6 7 8 9 10

7.621 9.810 9.948 11.135 11.216 12.482 12.872 13.650 14.676 14.738

7.410 9.726 9.764 10.896 11.114 12.353 12.781 13.368 14.485 14.766

7.548 10.764 11.113 11.328 12.862 13.300 14.168 15.369 16.205 17.137

7.428 9.829 9.858 10.960 11.178 12.367 12.835 13.433 14.440 14.743

applied in their approach to maintain C1 continuity between patches. In this work, we alternatively use the level sets to differentiate between inside and outside physical domains and the numerical integration is adopted only inside the physical domain, thus the process of partition NURBS patches and coupling patches are no longer required. The present results using 21  21 control points are listed in Tables 5–7 for completely free, simply supported and clamped boundary condition, respectively. The normalized frequencies ~ ¼ ½qhx2 a4 =D1=4 obtained by the proposed method are compared x with some available results derived from the IGA and Kirchhoff– Love theory with eight patches [33], MKI method [44], node-based

Table 8 Normalized frequencies of SSSS laminated square plate with complicated cutout for various orientations. Angle ply

Method

Mode 1

2

3

4

5

6

(0°, 0°, 0°)

This work IGA with 8-patches [30] EFG [47] MKI [47]

18.192 18.194 18.226 18.169

30.936 30.932 31.127 30.303

36.082 35.678 36.237 36.581

56.420 55.383 56.874 57.429

62.024 62.164 62.390 64.145

82.966 81.795 83.565 85.656

(15°, 15°, 15°)

This work IGA with 8-patches [30] EFG [47] MKI [47]

19.100 18.912 19.177 18.323

32.149 32.045 32.445 31.472

36.458 36.004 37.238 37.617

57.573 56.345 58.716 63.077

63.361 63.370 63.994 66.538

84.776 83.629 86.500 86.486

(30°, 30°, 30°)

This work IGA with 8-patches [30] EFG [47] MKI [47]

20.606 20.316 20.926 20.310

33.997 33.933 34.915 33.987

37.610 37.074 39.101 39.898

59.797 58.484 62.222 58.111

65.688 65.895 67.054 69.699

88.809 87.973 92.715 92.099

(45°, 45°, 45°)

This work IGA with 8-patches [30] EFG [47] MKI [47]

21.313 20.982 21.736 20.987

34.801 34.848 36.079 34.897

38.289 37.559 39.975 39.269

60.897 59.325 63.897 63.375

66.885 67.518 68.525 69.017

91.601 91.220 96.767 96.588

(0°, 90°, 0°)

This work IGA with 8-patches [30] EFG [47] MKI [47]

18.201 18.190 18.278 18.027

31.082 31.087 32.264 32.506

36.096 35.655 36.134 37.268

56.473 55.452 57.151 57.698

62.523 62.582 65.853 70.768

83.660 82.383 90.678 92.998

Table 9 Normalized frequencies of CCCC square plate with complicated cutout for various orientations. Angle ply

Mode 1

2

3

4

5

6

7

8

9

10

(0°, 0°, 0°) (15°, 15°, 15°) (30°, 30°, 30°) (45°, 45°, 45°) (0°, 90°, 0°)

44.269 45.084 45.736 45.998 44.428

71.432 71.802 73.395 76.075 72.232

81.620 81.245 79.268 76.253 81.787

97.102 96.875 96.545 96.393 97.180

101.593 101.969 101.865 101.819 101.823

119.213 119.492 121.494 125.184 119.684

135.113 133.960 130.752 126.208 135.070

147.394 144.010 140.870 139.979 147.578

162.798 165.086 168.343 170.452 163.539

171.948 172.305 172.999 173.072 172.669

160

S. Yin et al. / Composite Structures 127 (2015) 152–164

Table 10 Normalized frequencies of SCSC laminated square plate with complicated cutout for various orientations. Angle ply

(0°, 0°, 0°) (15°, 15°, 15°) (30°, 30°, 30°) (45°, 45°, 45°) (0°, 90°, 0°)

Mode 1

2

3

4

5

6

7

8

9

10

28.755 29.909 31.898 33.863 29.079

39.330 39.816 41.453 43.227 39.560

54.900 56.158 59.480 63.869 55.728

68.775 70.114 73.628 77.595 69.363

79.955 81.783 85.336 89.129 81.042

99.849 102.274 107.280 106.308 101.006

107.615 107.657 107.956 114.869 108.083

121.540 119.817 120.241 123.163 122.161

140.141 143.689 149.459 153.258 141.487

148.317 150.016 153.973 157.988 149.600

Table 11 First four mode shapes of three-layer symmetric laminated composite square plate with complicated cutout for different boundary conditions (B.C.).

smoothing RPIM (NS-RPIM) method [45] and EFG method [46]. It can be seen that the present results are in good agreement with the reference solutions for all cases of considered boundary conditions and, they close the results given by the IGA and Kirchhoff– Love theory with eight patches [33]. As compared to [33], fewer control points are required in this work and more importantly no any additional tasks, e.g., the bending strip method, are performed. As a result, the overall computational cost required in the present formation is less than the work proposed in [33]. 4.3. Composite plate with a complicated cutout We now turn to consider a three-layer symmetric laminated composite square plate with complicated cutout. The geometry of the plate and its dimensions are the same as those in Fig. 7. The material and geometrical parameters are considered as follows: ratio of elastic constants E1 =E2 ¼ 2:45 and G12 =E2 ¼ 0:48, Poisson’s ratio m12 ¼ 0:23, mass density q ¼ 8000 kg/m3 and thickness ~ ¼ ½qhx2 a4 =D0:1 1=2 h = 0.06 m. The frequency is normalized by x 3

in which a = 10 m and D0:1 ¼ E1 h =12ð1  v 12 v 21 Þ. The normalized first six frequencies of the three-ply laminated plate with SSSS boundary conditions are presented in Table 8 for various fiber orientations. The present solutions are also compared with results obtained by IGA and Kirchhoff with eight patches [30], EFG and MKI methods [47]. It clearly observed that the present results are in good agreement with the reference solutions for all considered orientations. Tables 9 and 10 tabulate the first ten normalized frequencies of CCCC and SCSC boundary conditions. The first four mode shapes of different boundary conditions of SSSS, SCSC and CCCC with fiber orientations (15°, 15°, 15°) are also shown in Table 11.

Fig. 8. Geometry parameters of a square plate with a complicated cutout.

The last benchmark example of eigenvalue analysis considers a four-layer square laminated plate with a complicated cutout as shown in Fig. 8. The material properties for the laminate plate

161

S. Yin et al. / Composite Structures 127 (2015) 152–164 Table 12 Normalized frequencies of CCCC laminated square plate with a hole of complicated cutout for various orientations. Angle ply

Method

Mode 1

2

3

4

5

6

(0°, 0°, 0°, 0°)

This work ABAQUS

29.623 29.573

32.573 32.565

47.009 46.767

51.321 51.237

67.057 66.712

69.474 69.463

(0°, 30°, 30°, 0°)

This work ABAQUS

31.182 31.155

34.800 34.731

48.267 48.143

56.120 56.052

71.022 70.838

71.337 71.182

(0°, 45°, 45°, 0°)

This work ABAQUS

32.968 32.909

37.695 37.826

49.259 49.174

61.773 61.554

72.562 72.558

75.529 75.309

(0°, 60°, 60°, 0°)

This work ABAQUS

34.955 34.731

41.280 41.265

50.044 49.862

67.271 67.056

74.287 73.933

80.032 79.779

(0°, 90°, 90°, 0°)

This work ABAQUS

37.103 37.023

45.655 45.661

50.623 50.416

72.562 72.187

76.279 76.053

84.402 84.189

Table 13 Normalized frequencies of SSSS laminated square plate with a hole of complicated cutout for various orientations. Angle ply

Method

Mode 1

2

3

4

5

6

(0°, 0°, 0°, 0°)

This work ABAQUS

8.339 8.322

12.287 12.276

16.710 16.609

21.999 21.939

35.121 34.731

39.859 39.890

(0°, 30°, 30°, 0°)

This work ABAQUS

10.174 10.144

14.781 14.752

20.531 20.426

28.370 28.301

38.557 38.514

44.633 44.704

(0°, 45°, 45°, 0°)

This work ABAQUS

10.671 10.660

15.515 15.509

22.275 22.180

32.709 32.599

39.189 38.858

47.338 47.111

(0°, 60°, 60°, 0°)

This work ABAQUS

10.188 10.179

15.106 15.096

22.794 22.696

35.712 35.419

38.433 38.170

48.580 48.486

(0°, 90°, 90°, 0°)

This work ABAQUS

8.613 8.600

13.620 13.627

22.285 22.149

36.355 36.167

37.610 37.460

48.411 48.336

FSM [8]

FSM [48]

FEM [49]

1.9837

2.000

2.090

Table 14 Comparisons of critical buckling coefficient for isotropic plate with square cutout. Method

This work Number of control points

k

99

16  16

21  21

28  28

33  33

2.0606

2.0177

1.999

1.9924

1.9886

(a)

(b)

Fig. 9. A square plate with a square cutout: geometric parameters (left) and mesh with 21  21 control points (right).

162

S. Yin et al. / Composite Structures 127 (2015) 152–164

Fig. 10. First eight mode shapes of a square plate with a square cutout.

are: ratio of elastic constants E1 =E2 ¼ 40 and G12 =E2 ¼ 0:6, Poisson’s ratio m12 ¼ 0:25, mass density q ¼ 8000 kg/m3 and thickness h = 0.01 m. The first six normalized frequencies ~ ¼ ½qhx2 a4 =D0:1 1=2 of various fiber orientations are listed in x Tables 12 and 13 for the two different boundary conditions, respectively. FEM results obtained by ABAQUS are also provided for comparison. The results of present method are in good agreement with the FEM using ABAQUS (see Table 14). 5. Buckling examples 5.1. Isotropic plate with a square cutout A square simply supported plate with central square cutout depicted in Fig. 9 is considered. The geometry and material properties are assumed as follows: a = 1 m, c/a = 0.5, thickness h = 0.01 m, Young’s modulus E ¼ 7:0  1010 N/m2, Poisson’s ratio m ¼ 0:3. The uniaxial in-plane compressive load is applied in the y-direction, and the critical buckling coefficient is normalized by 2

k ¼ 12ð1  v 2 Þða=hÞ N cr =p2 E. Five different meshes of 9  9, 16  16, 21  21, 28  28 and 33  33 control points are used. The comparison of the results from the present approaches with the corresponding results from the literature [8,48,49] shows the very good accuracy of the calculations. Fig. 10 plots the first eight buckling shape modes of this example. Fig. 11. Geometry of a square plate with a complicated cutout.

5.2. Isotropic plate with a complicated cutout The static buckling loads of a plate with complicated cutout are calculated using the present method. An in-plane compressive load is applied in the x-direction and the static buckling loads are represented via the factor of buckling load defined as k ¼ N cr a2 =p2 D. The geometry of the plate is shown in Fig. 11. The geometric parameters and material properties of the thin square plate are set as follows: length a = 10 m, thickness h = 0.06, Young’s modulus E ¼ 2:45  106 N/m2, Poisson’s ratio m ¼ 0:23. The effects of different boundary conditions on the buckling load factor are investigated. Table 15 presents the buckling load factors obtained by the developed method employing cubic NURBS basis function with a 21  21 control mesh. Once again, a good agreement between the proposed approach and reference solutions [50] is found.

Table 15 Buckling factor k for the plate with a complicated cutout. Boundary condition

EFG [50]

This work

SSSS CCCC SCSC SCCS

2.64 16.66 5.43 3.97

2.60 17.41 5.43 3.95

5.3. Composite plate with a complicated cutout Last example considers a three-layer symmetric laminated composite plate with a complicated cutout as depicted in Fig. 11. The geometry of the plate and its dimensions are the same as those

163

S. Yin et al. / Composite Structures 127 (2015) 152–164 Table 16 Buckling loads of a SSSS plate with complicated cutout for various orientations. Angle ply

Mode

(0°, 0°, 0°) (15°, 15°, 15°) (30°, 30°, 30°) (45°, 45°, 45°) (0°, 90°, 0°)

1

2

3

4

5

6

7

8

9

10

1.360 1.479 1.715 1.827 1.359

2.319 2.402 2.563 2.603 2.317

3.448 3.706 4.221 4.495 3.474

5.689 5.963 6.507 6.766 5.710

8.974 8.886 8.508 7.790 8.959

12.166 12.034 11.416 10.338 12.131

12.556 12.307 11.604 10.632 12.416

15.007 14.450 13.185 11.760 14.808

17.986 18.357 19.041 18.942 18.003

18.665 18.733 19.064 18.987 18.642

Table 17 Buckling loads of a CCCC plate with complicated cutout for various orientations. Angle ply

Mode

(0°, 0°, 0°) (15°, 15°, 15°) (30°, 30°, 30°) (45°, 45°, 45°) (0°, 90°, 0°)

1

2

3

4

5

6

7

8

9

10

11.240 11.053 10.798 10.513 11.294

11.436 11.184 10.925 10.610 11.512

11.674 11.791 11.725 11.551 11.688

11.894 11.951 11.895 11.698 11.929

22.627 21.223 18.311 15.649 22.311

23.065 21.534 18.537 15.818 22.711

24.131 23.162 20.533 17.731 23.840

24.680 23.553 20.783 17.894 24.337

34.625 33.631 31.132 28.078 34.427

35.391 34.189 31.402 28.207 35.123

Table 18 Buckling loads of a SCSC plate with complicated cutout for various orientations. Angle ply

Mode 1

2

3

4

5

6

7

8

9

10

(0°, 0°, 0°) (15°, 15°, 15°) (30°, 30°, 30°) (45°, 45°, 45°) (0°, 90°, 0°)

2.861 2.971 3.234 3.476 2.902

2.925 3.032 3.279 3.497 2.960

7.597 7.923 8.755 9.684 7.759

7.786 8.113 8.914 9.763 7.928

13.828 13.572 12.991 12.345 13.874

14.359 14.076 13.360 12.512 14.365

15.564 15.178 14.336 13.460 15.466

16.081 15.628 14.625 13.583 15.941

22.004 22.320 22.875 22.440 22.178

22.204 22.471 22.968 22.543 22.362

in Fig. 11. The material and geometrical parameters are considered as follows: ratio of elastic constants E1 =E2 ¼ 2:45 and G12 =E2 ¼ 0:48, Poisson’s ratio m12 ¼ 0:23 and thickness h = 0.06 m. The buckling load coefficient is normalized by k ¼ N cr a2 =p2 D0:1 in 3

which D0:1 ¼ E1 h =12ð1  v 12 v 21 Þ. The effects of different boundary conditions and fiber orientations on the buckling load coefficient are considered and the obtained numerical results are presented in Tables 16–18. For the SSSS boundary condition, the maximum buckling load coefficient of the first, second, third and fourth modes occur at h ¼ 45 , the fifth, sixth, seventh and eighth modes occur at h ¼ 0 , the ninth and tenth modes occur at h ¼ 30 . For the CCCC boundary condition, maximum buckling load coefficient of the first and second, modes occur at (0°, 90°, 0°), the third and fourth modes occur at h ¼ 15 , the others modes occur at h ¼ 0 . For the SCSC boundary condition, the maximum buckling load coefficient of the first, second, third and fourth modes occur at h ¼ 45 , the fifth and sixth modes occur at (0°, 90°, 0°), the seventh and eighth modes occur at h ¼ 0 , the ninth and tenth modes occur at h ¼ 30 . It also can be seen that the buckling loads with clamped boundaries are generally higher than those with simply supported boundaries. 6. Conclusions We develop an effective and accurate NURBS-based IGA using the classical plate theory for the analysis of free vibration and buckling behaviors of thin laminated composite structures with cutouts. The C1 consistency requirement of classical plate theory is straightforward to construct due to the higher-order continuity of NURBS basis function. However, the tensor product of the NURBS basis functions induces the difficulty in treating the

trimmed objects like internal cutouts. In this work, the level sets are used to describe the cutouts and the numerical integration is used only inside the physical domain. As a result, the trimmed NURBS surface is thus no longer required to describe the geometrical structure with cutouts. The accuracy and efficiency of the proposed method are illustrated through numerical examples with different complicated geometries and boundary conditions. Compared results have shown that the proposed NURBS-based IGA analysis is effective and well suited for dynamic and buckling analysis of thin laminates with cutouts. Acknowledgments This work was supported by the scholarship from China Scholarship Council (CSC), Jiangsu Province Graduate Students Research and Innovation Plan (Grant No. CXZZ13_0235), the National Natural Science Foundation of China (Grant No. 51179063), the National Sci-Tech Support Plan of China (Grant No. 2015BAB07B10) and the Grant-in-Aid for Scientific Research (No. 26-04055) – Japan Society for the Promotion of Sciences (JSPS, ID No. P14055). The financial supports are gratefully acknowledged. References [1] Chai GB. Free vibration of laminated composite plates with a central circular hole. Compos Struct 1997;35:357–68. [2] Sivakumar K, Iyengar NGR, Deb K. Free vibration of laminated composite plates with cutouts. J Sound Vib 1999;221(3):443–70. [3] Aydin Komur M, Sen Faruk, Atas Akin, Arslan Nurettin. Buckling analysis of laminated composite plates with an elliptical/circular cutout using fem. Adv Eng Softw 2010;41:161–4. [4] Kumar D, Singh SB. Effects of boundary conditions on buckling and postbuckling responses of composite laminate with various shaped cutouts. Compos Struct 2010;92:769–79.

164

S. Yin et al. / Composite Structures 127 (2015) 152–164

[5] Chen CC, Kitiporanchai S, Lim CW, Liew KM. Free vibration of symmetrically laminated thick perforated plates. J Sound Vib 2000;230:111–32. [6] Liu GR, Zhao X, Dai KY, Zhong ZH, Li GY, Han X. Static and free vibration analysis of laminated composite plates using the conforming radial point interpolation method. Compos Sci Technol 2008;68:354–66. [7] Eccher G, Rasmussen KJR, Zandonini R. Elastic buckling analysis of perforated thin-walled structures by the isoparametric spline finite strip method. ThinWalled Struct 2008;46:165–91. [8] Ovesy HR, Fazilati J. Buckling and free vibration finite strip analysis of composite plates with cutout based on two different modeling approaches. Compos Struct 2012;94(3):1250–8. [9] Natarajan S, Deogekar PS, Manickam G, Belouettar S. Hygrothermal effects on the free vibration and buckling of laminated composites with cutouts. Compos Struct 2014;108:848–55. [10] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng 2005;194(39-41):4135–95. [11] Yu TT, Yin SH, Bui QT, Hirose S. A simple FSDT-based isogeometric analysis for geometrically nonlinear analysis of functionally graded plates. Finite Elem Anal Des 2015;96:1–10. [12] Bazilevs Y, Beirao Da Veiga L, Cottrell JA, Hughes TJR, Sangalli G. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes. Math Models Methods Appl Sci 2006;16(07):1031–90. [13] Evans JA, Bazilevs Y, Babuška I, Hughes TJR. n-Widths, sup–infs, and optimality ratios for the k-version of the isogeometric finite element method. Comput Methods Appl Mech Eng 2009;198(21-26):1726–41. [14] Hughes TJR, Reali A, Sangalli G. Efficient quadrature for NURBS-based isogeometric analysis. Comput Methods Appl Mech Eng 2010;199(58):301–13. [15] Auricchio F, Calabrò F, Hughes TJR, Reali A, Sangalli G. A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis. Comput Methods Appl Mech Eng 2012;249-252:15–27. [16] Bazilevs Y, Calo VM, Cottrell JA, Evans JA, Hughes TJR, Lipton S, et al. Isogeometric analysis using T-splines. Comput Methods Appl Mech Eng 2010;199(5-8):229–63. [17] Valizadeh N, Natarajan S, Gonzalez-Estrada OA, Rabczuk T, Bui TQ, Bordas SPA. NURBS-based finite element analysis of functionally graded plates: static bending, vibration, buckling and flutter. Compos Struct 2013;99:309–26. [18] Valizadeh N, Bui QT, Vu TV, Thai HT, Nguyen NM. Isogeometric simulation for buckling, free and forced vibration of orthotropic plates. Int J Appl Mech 2013;05:1350017. [19] Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. Isogeometric shell analysis: the Reissner–Mindlin shell. Comput Methods Appl Mech Eng 2010;199(58):276–89. [20] Yin SH, Hale JS, Yu TT, Bui TQ, Bordas SPA. Isogeometric locking-free plate element: a simple first order shear deformation theory for functionally graded plates. Compos Struct 2014;118:121–38. [21] Bazilevs Y, Calo VM, Cottrell JA, Hughes TJR, Reali A, Scovazzi G. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput Methods Appl Mech Eng 2007;197(14):173–201. [22] Bazilevs Y, Calo VM, Hughes TJR, Zhang Y. Isogeometric fluid–structure interaction: theory, algorithms, and computations. Comput Mech 2008;43(1):3–37. [23] Verhoosel CV, Scott MA, Hughes TJR, de Borst R. An isogeometric analysis approach to gradient damage models. Int J Numer Methods Eng 2011;86(1):115–34. [24] Lu J. Isogeometric contact analysis: geometric basis and formulation for frictionless contact. Comput Methods Appl Mech Eng 2011;200:726–41. [25] Nguyen MN, Bui QT, Yu TT, Hirose S. Isogeometric analysis for unsaturated flow problems. Comput Geotech 2014;62:257–67. [26] Wall WA, Frenzel MA, Cyron C. Isogeometric structural shape optimization. Comput Methods Appl Mech Eng 2008;197(33–40):2976–88.

[27] Kiendl J, Bletzinger KU, Linhard J, Wüchner R. Isogeometric shell analysis with Kirchhoff–Love elements. Comput Methods Appl Mech Eng 2009;198(4952):3902–14. [28] Kiendl J, Bazilevs Y, Hsu MC, Wüchner R, Bletzinger KU. The bending strip method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches. Comput Methods Appl Mech Eng 2010;199(37– 40):2403–16. [29] Benson DJ, Bazilevs Y, Hsu MC, Hughes TJR. A large deformation, rotation-free, isogeometric shell. Comput Methods Appl Mech Eng 2011;200(13– 16):1367–78. [30] Shojaee S, Valizadeh N, Izadpanah E, Bui TQ, Vu TV. Free vibration and buckling analysis of laminated composite plates using the NURBS-based isogeometric finite element method. Compos Struct 2012;94(5):1677–93. [31] Yin SH, Yu TT, Liu P. Free vibration analyses of FGM thin plates by isogeometric analysis based on classical plate theory and physical neutral surface. Adv Mech Eng 2013. Article ID 634584, 10p. [32] Lu J, Zheng C. Dynamic cloth simulation by isogeometric analysis. Comput Methods Appl Mech Eng 2014;268:475–93. [33] Shojaee S, Izadpanah E, Valizadeh N, Kiendl J. Free vibration analysis of thin plates by using a NURBS-based isogeometric approach. Finite Elem Anal Des 2012;61:23–34. [34] Schmidt R, Wüchner R, Bletzinger KU. Isogeometric analysis of trimmed NURBS geometries. Comput Methods Appl Mech Eng 2012;241-244:93–111. [35] Ghorashi SS, Valizadeh N, Mohammadi S. Extended isogeometric analysis for simulation of stationary and propagating cracks. Int J Numer Methods Eng 2011;89(9):1069–101. [36] Rank E, Ruess M, Kollmannsberger S, Schillinger D, Düster A. Geometric modeling, isogeometric analysis and the finite cell method. Comput Methods Appl Mech Eng 2012;249–252:104–15. [37] Schillinger D, Dede’ L, Scott MA, Evans JA, Borden MJ, Rank E, et al. An isogeometric design-through-analysis method based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces. Comput Methods Appl Mech Eng 2012;249–252:116–50. [38] Piegl L, Tiller W. The NURBS book. Berlin Heidelberg: Springer; 1995. [39] Bhardwaj G, Singh IV, Mishra BK, Bui QT. Numerical simulation of functionally graded cracked plates using NUBRS based XIGA under different loads and boundary conditions. Compos Struct; 2015, in press, http://dx.doi.org/10. 1016/j.compstruct.2015.02.066. [40] Bui QT, Zhang Ch. Analysis of generalized dynamic intensity factors of cracked magnetoelectroelastic solids by X-FEM. Finite Elem Anal Des 2013;69:19–36. [41] Bui QT, Zhang Ch. Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading. Comput Mater Sci 2012;62:243–57. [42] Nagy AP, Benson DJ. On the numerical integration of trimmed isogeometric elements. Comput Methods Appl Mech Eng 2015;284:165–85. [43] Huang M, Sakiyama T. Free vibration analysis of rectangular plates with variously-shaped holes. J Sound Vib 1999;226:769–86. [44] Bui TQ, Nguyen MN. A moving Kriging interpolation-based meshfree method for free vibration analysis of Kirchhoff plates. Comput Struct 2011;89(3– 4):380–94. [45] Cui XY, Liu GR, Li GY, Zhang GY. A thin plate formulation without rotation DOFs based on the radial point interpolation method and triangular cells. Int J Numer Methods Eng 2011;85(8):958–86. [46] Liu GR, Chen XL. A mesh-free method for static and free vibration analyses of thin plates of complicated shape. J Sound Vib 2001;241(5):839–55. [47] Bui QT, Nguyen NM, Zhang Ch. An efficient meshfree method for vibration analysis of laminated composite plates. Comput Mech 2011;48:175–93. [48] Cheung YK, Kong J. Linear elastic stability analysis of shear-deformable plates using a modified spline finite strip method. Compos Struct 1993;41(2):189–92. [49] Tham LG, Chan AHC, Cheung YK. Free vibration and buckling analysis of plates by the negative stiffness method. Compos Struct 1986;22(4):687–92. [50] Liu GR. Meshfree methods: moving beyond the finite element method. 2nd ed. CRC Press; 2010.