Transient analysis of laminated composite plates using nurbs-Based isogeometric analysis

Volume 36 Number 4 4 2014 Vietnam Journal of Mechanics, VAST, Vol. 36, No. 4 (2014), pp. 267 – 281 TRANSIENT ANALYSIS OF LAMINATED COMPOSITE PLATES USING NURBS-BASED ISOGEOMETRIC ANALYSIS Lieu B. Nguyen1, Chien H. Thai2, Ngon T. Dang1, H. Nguyen-Xuan3,∗ 1Ho Chi Minh City University of Technology and Education, Vietnam 2Ton Duc Thang University, Ho Chi Minh City, Vietnam 3Vietnamese-German University, Ho Chi Minh City, Vietnam ∗E-mail: hung.nx@vgu.edu.vn Received July 04, 2014 Abstract. We further study isogeometric approach for response analysis of laminated composite plates using the higher-order shear deformation theory. The present theory is derived from the classical plate theory (CPT) and the shear stress free surface conditions are naturally satisfied. Therefore, shear correction factors are not required. Galerkin weak form of response analysis model for laminated composite plates is used to obtain the discrete system of equations. It can be solved by isogeometric approach based on the non-uniform rational B-splines (NURBS) basic functions. Some numerical examples of the laminated composite plates under various dynamic loads, fiber orientations and lay-up numbers are provided. The accuracy and reliability of the proposed method is verified by comparing with analytical solutions, numerical solutions and results from Ansys software. Keywords: Transient analysis, laminated composite plate, isogeometric analysis, NURBS, Newmark integration. 1. INTRODUCTION The transient response of laminated composite plates has received much attention from designers due to increasing applications of composite in high performance aircraft, vehicles and vessels. Whether they are used in civil, marine or aerospace, most structures are subjected to dynamic loads during their operation. Therefore, there exists a need for assessing the natural frequency and transient response of structures. Many numerical methods have been developed to compute, analyze and simulate the response as well as dynamic characteristics of laminated composite plates. Out of these methods, the finite element method (FEM) has become the universally applicable technique for solving boundary and initial value problems. In the past years, Reismann [1], Reismann and Lee [2] have analyzed simply supported rectangular isotropic plates, which are subjected to suddenly a uniformly distributed load over a square area on the plate. The transient finite element analysis of isotropic plate was also carried out by 268 Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan Rock and Hinton [3] for thick and thin plates. Akay [4] determined the large deflection transient response of isotropic plates using a mixed FEM. For composite plates, Reddy [5] presented finite element results for the transient analysis of layered composite plates based on the first-order shear deformation theory (FSDT). Mallikarjuna and Kant [6] presented an isoparametric finite element formulation based on a higher-order displacement model for dynamic analysis of multi-layer symmetric composite plate. Wang and his co-workers developed the strip element method (SEM) for static bending analysis of orthotropic plates. Then, Wang et al. [7] extended the SEM to analyze dynamic response of symmetric laminated plates. Although FEM is an extremely versatile and powerful technique, it has certain dis- advantages. Recently, Hughes and his co-workers have proposed a robustly computational isogeometric analysis [8]. Following this approach, the CAD-shape functions, commonly the non-uniform rational B-splines (NURBS) are substituted for the Lagrange polynomial based shape functions in the CAE. The computational cost is decreased significantly as the meshes are generated within the CAD. IGA gives higher accurate results because of the smoothness and the higher-order continuity between elements [9, 10]. In this paper, a higher-order displacement field in which the in-plane displacement is expressed as cubic functions of the thickness coordinate with constant transverse displace- ment across the thickness is used. The finite element formulation based on the higher-order shear deformation theory (HSDT) requires elements with at least C1-inter-element con- tinuity. It is difficult to achieve such elements for free-form geometries when using the standard Lagrangian polynomials as basis functions. Fortunately, IGA can be easily ob- tained because NURBS basis functions are Cp−1 continuous. The governing equations of the laminated composite plates are transformed into a standard weak-form, which is then discretized into the system of time-dependent equations to be solved by the unconditionally stable Newmark time integration scheme. Several numerical examples with many different models are provided to illustrate the effectiveness and reliability of the present method in comparison with other results from the literature. The paper is outlined as follows. Next section introduces the HSDT for laminated composite plates. In section 3, the numerical formulation relied on the HSDT and IGA is described. The numerical results and discussions are provided in section 4. Finally, in section 5, concluding remarks are presented with the brief discussion of the numerical results obtained by the developed methodology. 2. THE HIGHER-ORDER SHEAR DEFORMATION THEORY FOR PLATES Let Ω be the domain in R2 occupied by the mid-plane of the plate and u0, v0, w and β = (βx;βy)T denote the displacement components in the x; y; z directions and the rotations in the x−z and y−z planes (or the-y and the-x axes), respectively. Fig. 1 shows the geometry of a plate and the coordinate system. A generalized displacement field of an arbitrary point in the plate based on higher-order shear deformation theory derived from the classical plate theory is defined as follows [9] Transient analysis of laminated composite plates using NURBS-based isogeometric analysis 269 u (x, y, z, t) = u0 (x, y, t)− z ∂w (x, y, t) ∂x + f (z)βx (x, y, t) v (x, y, z, t) = v0 (x, y, t)− z ∂w (x, y, t) ∂y + f (z)βy (x, y, t), (−h 2 ≤ z ≤ h 2 ) w (x, y, z, t) = w (x, y, t) (1) In this paper we exploit the third-order shear deformation theory (TSDT) of Reddy [11] and the distribution function is written as f (z) = z − 4z3/3h2. 2 elements for free-form geometries when using the standard Lagrangian polynomials as basis functions. Fortunately, IGA can be easily obtained because NURBS basis functions are Cp-1 continuous. The governing equations of the laminated composite plates are transformed into a standard weak-form, which is then discretized into the system of time-dependent equations to be solved by the unconditionally stable Newmark time integration scheme. Several numerical examples with many different models are provided to illustrate the effectiveness and reliability of the present method in comparison with other results from the literature. The paper is outlined as follows. Next section introduces the HSDT for laminated composite plates. In section 3, the numerical formulation relied on the HSDT and IGA is described. The numerical results and discussions are provided in section 4. Finally, in section 5, concluding remarks are presented with the brief discussion of the numerical results obtained by the developed methodology. 2. THE HIGHER-ORDER SHEAR DEFORMATION THEORY FOR PLATES Let  be the domain in R2 occupied by the mid-plane of the plate and u0, v0, w and  = (x ;y) T denote the displacement components in the x; y; z directions and the rotations in the x-z and y-z planes (or the-y and the-x axes), respectively. Fig. 1 shows the geometry of a plate and the coordinate system. A generalized displacement field of an arbitrary point in the plate based on higher-order shear deformation theory derived from the classical plate theory is defined as follows [9]:                         0 0 , , , , , , , , , , , , , , , , , , , , , , , x y w x y t u x y z t u x y t z f z x y t x w x y t v x y z t v x y t z f z x y t y w x y z t w x y t              ; 2 2 h h z        (1) In this paper, we exploit the third-order shear deformation theory (TSDT) of Reddy [11] and the distribution function is written as   3 24 3  /f z z z h . Fig. 1. Plate model and coordinate system. The relationship between strains and displacements is described by, 0 1 2[ ] ( ) T p xx yy xy z f z         ( ) T xz yz sf z     γ ε (2) where Fig. 1. Plate model and coordinate system Th relationship between str ins and di placements is described by εp = [εxx εyy γxy] T = ε0 + zε1 + f(z)ε2, γ = [γxz γyz] T = f ′(z)εs (2) where ε0 =  ∂u0 ∂x ∂v0 ∂y ∂v0 ∂x + ∂u0 ∂y  , ε1 =  −∂ 2w ∂x2 −∂ 2w ∂y2 −2 ∂ 2w ∂x∂y  , ε2 =  ∂βx ∂x ∂βy ∂y ∂βy ∂x + ∂βx ∂y  , εs = [ βx βy ] (3) Neglecting σz for each orthotropic layer, the constitutive equation of an orthotropic lamina in the local coordinate system is derived from Hooke’s law for a plane stress con- dition as  σk1 σk2 τk12 τk13 τk23  =  Q11 Q12 0 0 0 Q12 Q22 0 0 0 0 0 Q33 0 0 0 0 0 Q55 0 0 0 0 0 Q44  k εk1 εk2 γk12 γk13 γk23  , (4) 270 Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan where subscripts 1 and 2 are the directions of the fiber and in-plane normal to fiber, respectively, subscript 3 indicates the direction normal to the plate, and the reduced stiffness components, Qkij are given by Qk11 = Ek1 1− νk12νk21 , Qk12 = νk12E k 2 1− νk12νk21 , Qk22 = Ek2 1− νk12νk21 , Qk33 = G k 12, Q k 55 = G k 13, Q k 44 = G k 23, (5) in which E1, E2, G12, G23, G13 and ν12 are independent material properties for each layer. The laminate is usually made of several orthotropic layers. Each layer must be transformed into the laminate coordinate system (x, y, z). The stress-strain relationship is given as  σkxx σkyy τkxy τkxz τkyz  =  Q¯11 Q¯12 Q¯16 0 0 Q¯12 Q¯22 Q¯26 0 0 Q¯61 Q¯62 Q¯33 0 0 0 0 0 Q¯55 Q¯54 0 0 0 Q¯45 Q¯44  k εkxx εkyy γkxy γkxz γkyz  , (6) where Q¯kij is the transformed material constant matrix (see [12] for more details). From Hooke’s law and the linear strains given by Eq. (2), the stress is computed by σ = [ σp τ ] = [ D∗ 0 0 Ds ] [ εp γ ] , (7) where σp and τ are the in-plane stress component and shear stress, respectively, and D∗ is material constant matrices given in the form as D∗ =  A B EB D F E F H  , (8) where Aij , Bij , Dij , Eij , Fij , Hij = ∫ h/2 −h/2 (1, z, z2, f(z), zf(z), f2(z))Q¯ijdz, i, j = 1, 2, 6, Dsij = ∫ h/2 −h/2 [ (f ′(z))2 ] Q¯ijdz, i, j = 4, 5. (9) For forced vibration analysis of the plates, a weak form can be derived from the following undamped dynamic equilibrium equation as∫ Ω δεTpD ∗εpdΩ + ∫ Ω δγTDsγdΩ + ∫ Ω δu˜Tm¨˜udΩ = ∫ Ω δwq(x, y, t)dΩ, (10) where q(x, y, t) is the transverse loading per unit area and the function depending on time and space; the mass matrix m is calculated according to the consistent form given by m =  I1 I2 I4I2 I3 I5 I4 I5 I6  , (I1, I2, I3, I4, I5, I6) = h/2∫ −h/2 ρ ( 1, z, z2, f(z), zf(z), f2(z) ) dz, (11) Transient analysis of laminated composite plates using NURBS-based isogeometric analysis 271 in which u˜ = [ u1 u2 u3 ]T with u1 =  u0v0 w  , u2 =  −w,x−w,y 0  , u3 =  βxβy 0  , (12) where ρ is the mass density. 3. THE LAMINATED COMPOSITE PLATE FORMULATION BASED ON NURBS BASIS FUNCTIONS 3.1. Introduction to NURBS basis functions [9] Given a knot vector Ξ = {ξ1, ξ2, . . . , ξn+p+1}, the associated B-spline basis functions are defined recursively starting with the zeroth order basis function (p = 0) as Ni,0 (ξ) = { 1 if ξi ≤ ξ < ξi+1 0 otherwise } , (13) and for a polynomial order p ≥ 1 Ni,p (ξ) = ξ − ξi ξi+p − ξiNi,p−1 (ξ) + ξi+p+1 − ξ ξi+p+1 − ξi+1Ni+1,p−1 (ξ) . (14) A knot vector Ξ is defined as a sequence of knot value ξi ∈ R, i = 1, . . . , n + p. If the first and the last knots are repeated p+ 1 times, the knot vector is called open knot. By the tensor product of basis functions in two parametric dimensions ξ and η with two knot vectors Ξ = {ξ1, ξ2, . . . , ξn+p+1} and H = {η1, η2, . . . , ηm+q+1}, the two- dimensional B-spline basis functions are obtained as, NA (ξ, η) = Ni,p (ξ)Mj,q (η). Fig. 2 illustrates a bivariate cubic B-spline basic function. 5 Fig. 2. A bivariate cubic B-spline basis function with knot vectors  0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1 Ξ Η 3.2. A higher order plate formulation based on NURBS approximation Using the NURBS basis functions defined above, both the description of the geometry (or the physical point) and the displacement field u of the plate are approximated as,     1 , , ; m n h A A A R      x P     1 , , m n h A A A R      u q (15) where n×m is the number basis functions,  T x yx is the physical coordinate vector. In Eq. (15),  ,AR   is rational basic functions, AP is the control points and 0 0 T A A A A xA yAu v w     q is the vector of nodal degrees of freedom associated with the control point A. Substituting Eq. (15) into Eq. (3), the in-plane and shear strains can be rewritten as: 1 2 1 m n T T m b b s p A A A A A A          B B B B q  (16) in which , , 1 , , , , , 0 0 0 0 0 0 0 0 0 0 0 0 , 0 0 0 0 0 0 0 0 0 2 0 0 A x A xx m b A A y A A yy A y A x A xy R R R R R R R                      B B , 2 , , , 0 0 0 0 0 0 0 0 0 0 0 0 , 0 0 0 0 0 0 0 A x Ab s A A y A A A y A x R R R R R R                B B (17) For forced vibration analysis of the plates, undamped dynamic discrete equation can be written from Eq. (10) as, (t)Mq + Kq = F (18) where the global stiffness matrix K is given by Fig. 2. A bivariate cubic B-spline basis function with knot vectors Ξ = H = {0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1} 272 Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan To exactly represent some curved geometries (e.g. circles, cylinders, spheres, etc.) the non-uniform rational B-splines (NURBS) functions are used. Being different from B- spline, each control point of NURBS has additional value called an individual weight ζA [8]. Thus, the NURBS functions can be expressed as RA (ξ, η) = NAζA/ m×n∑ A=1 NA (ξ, η) ζA. It is clear that B-spline function is only the special case of the NURBS function when the individual weight of control point is constant. 3.2. A higher-order plate formulation based on NURBS approximation Using the NURBS basis functions defined above, both the description of the geom- etry (or the physical point) and the displacement field u of the plate are approximated as xh (ξ, η) = m×n∑ A=1 RA (ξ, η)PA; uh (ξ, η) = m×n∑ A=1 RA (ξ, η)qA, (15) where n×m is the number basis functions, xT = (x y) is the physical coordinate vector. In Eq. (15), RA (ξ, η) is rational basic functions, PA is the control points and qA = [ u0A v0A wA βxA βyA ]T is the vector of nodal degrees of freedom associ- ated with the control point A. Substituting Eq. (15) into Eq. (3), the in-plane and shear strains can be rewritten as [εp γ] T = m×n∑ A=1 [ BmA B b1 A B b2 A B s A ]T qA, (16) in which BmA =  RA,x 0 0 0 00 RA,y 0 0 0 RA,y RA,x 0 0 0  , Bb1A =  0 0 −RA,xx 0 00 0 −RA,yy 0 0 0 0 −2RA,xy 0 0  Bb2A =  0 0 0 RA,x 00 0 0 0 RA,y 0 0 0 RA,y RA,x  , BsA = [ 0 0 0 RA 00 0 0 0 RA ] . (17) For forced vibration analysis of the plates, undamped dynamic discrete equation can be written from Eq. (10) as Mq¨ +Kq = F(t), (18) where the global stiffness matrix K is given by K = ∫ Ω   Bm Bb1 Bb2  T  A B EB D F E F H  Bm Bb1 Bb2 + (Bs)TDsBs dΩ. (19) The distributed transverse force in the z direction one has the following expression F(t) = ∫ Ω Rq(x, y, t)dΩ, (20) Transient analysis of laminated composite plates using NURBS-based isogeometric analysis 273 where R = [ 0 0 RA 0 0 ] . (21) The global mass matrix M is given as M = ∫ Ω   N1N2 N3 T  I1 I2 I4I2 I3 I5 I4 I5 I6   N1N2 N3   dΩ, (22) where N1 = RA 0 0 0 00 RA 0 0 0 0 0 RA 0 0 ;N2 = 0 0 −RA,x 0 00 0 −RA,y 0 0 0 0 0 0 0 ;N3 = 0 0 0 RA 00 0 0 0 RA 0 0 0 0 0 . (23) It should be noted that for forced vibration analysis, the approximate function is done with both space and time. For the displacements and accelerations at time t + ∆t, Eq. (18) should be considered at time t+ ∆t as follows Mq¨t+∆t + Kqt+∆t = Ft+∆t(t). (24) To solve this second order time dependent problem, several methods have been proposed such as, Wilson, Newmark, Houbolt, Crank-Nicholson, etc. In this paper, Eq. (18) is solved by the Newmark time integration method. The Newmark method is an implicit method. This method assumes that the acceleration varies linearly within the interval (t, t+ ∆t). The formulation of the Newmark method is [13][ M + αK(∆t)2 ] q¨1 = F1 − [Kq0 + K∆tq˙0 + ( 1 2 − α)Kq¨0(∆t)2], (25) q˙1 = q˙0 + (1− δ)q¨0∆t+ δq¨1∆t, (26) q1 = q0 + q˙0∆t+ ( 1 2 − α)q¨0(∆t)2 + αq¨1(∆t)2. (27) The parameters α and δ are constants whose values depend on the finite difference scheme used in the calculations. Two well-known and commonly used cases are average acceleration method (α = 1/4 and δ = 1/2) and linear acceleration method (α = 1/6 and δ = 1/2). Here we used the average acceleration method, which is unconditionally stable if δ ≥ 0.5 and α ≥ 1/4(δ + 0.5)2. 4. NUMERICAL EXAMPLES 4.1. A study of the convergence Free vibration analysis of the laminated composite plates is investigated correspond- ing to right hand side of Eq. (18) equal to zero. Let us consider a four-layer [00/900/900/00] square plate with simply supported boundary condition. The effects of the length to thick- ness a/h and elastic modulus ratios E1/E2 are studied. To show the convergence of the present approach, the length to thickness a/h = 5 and elastic modulus ratios E1/E2 = 40 are used. As shown in Tab. 1, the normalized frequencies are computed using meshes of 274 Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan 9 × 9, 13 × 13 and 17 × 17. It can be observed that the differences of normalized fre- quencies between meshes of 9 × 9 and 13 × 13 are not significant and between meshes of 13 × 13 and 17 × 17 are identical. Hence, for a comparison with other methods, a mesh of 13× 13 cubic elements can be chosen. The first normalized frequency derived from the present approach is listed in Tab. 2 corresponding to various modulus ratios and a/h = 5. The obtained results are compared with analytical solutions based on the HSDT [14, 15] the moving least squares differential quadrature method (DQM) [16] based on the FSDT, the meshfree method using multiquadric radial basis functions (RBFs) [17] and wavelets functions [18] based on the FSDT. A good agreement is found for the present method in comparison with other ones. It is also seen that the present results match very well with the exact solutions [14, 15]. The influence of the length to thickness ratios is also consid- ered as displayed in Tab. 3. The obtained results are compared with those of Zhen and Wanji [19] based on a global-local higher-order theory (GLHOT), Matsunaga [20] based on a global-local higher-order theory. Again, a good agreement with other published solutions is obtained. Table 1. The convergence of non-dimensional frequency parameter $ = ( ωa2/h ) (ρ/E2) 1/2 of a four layer [00/900/900/00] simply supported laminated square plate (a/h = 5) Method Meshes 9× 9 13× 13 17× 17 IGA (present) 10.7876 10.7873 10.7873 Table 2. A non-dimensional frequency parameter $ = ( ωa2/h ) (ρ/E2) 1/2 of a [00/900/900/00] simply supported laminated square plate (a/h = 5) Method E1/E2 10 20 30 40 RBFs-FSDT [17] 8.2526 9.4974 10.2308 10.7329 Wavelets-FSDT [18] 8.2794 9.5375 10.2889 10.8117 DQM-FSDT [16] 8.2924 9.5613 10.3200 10.8490 Exact-HSDT [14,15] 8.2982 9.5671 10.3260 10.8540 IGA (present) 8.2718 9.5263 10.2719 10.7873 Table 3. A non-dimensional frequency parameter $ = ( ωa2/h ) (ρ/E2) 1/2 of a [00/900/900/00] simply supported laminated square plate (E1/E2 = 40) Methods a/h 4 5 10 20 25 50 100 Zhen et al. [19] 9.2406 10.7294 15.1658 17.8035 18.2404 18.9022 19.1566 Matsunaga [20] 9.1988 10.6876 15.0721 17.6369 18.0557 18.6702 18.8352 IGA (present) 9.3235 10.7873 15.1073 17.6466 18.0620 18.6718 18.8356 Transient analysis of laminated composite plates using NURBS-based isogeometric analysis 275 4.2. Transient analysis In order to demonstrate the accuracy and effectiveness of the present method for transient analysis of laminated composite plates, four numerical examples with different transient loadings are studied. The obtained results are compared with other numerical or analytical solutions available in the literature or commercial software. For first three examples, cubic order NURBS basis function with 13×13 elements is used. All layers of the laminated plates are assumed to have the same thicknesses and material properties. The time step ∆t = 0.1 ms is chosen for Sections 4.2.1 and 4.2.2. 4.2.1. A three-layer square plate [00/900/00] First, a fully simply supported three-layer square laminated plate sorted as [00/900/00] is considered. Material I is used, shown in Tab. 4. This example was also studied by Wang et al. [7] using the trip element method (SEM), which is chosen here to demonstrate the accuracy of the IGA in dynamic analysis of plates under different tran- sient loads including step, triangular, sine and explosive blast loads. The length and the thickness of square plate are assumed to be a = 20h and h = 0.0381 m, respectively. The plate is subjected to a transverse load which is sinusoidally distributed in spatial domain and varies with time as q(x, y, t) = q0 sin( pix a ) sin( piy b )F (t), (28) in which F (t) =  { 1 0 ≤ t ≤ t1 0 t > t1 } Step loading{ 1− t/t1 0 ≤ t ≤ t1 0 t > t1 } Triangular loading{ sin(pit/t1) 0 ≤ t ≤ t1 0 t > t1 } Sine loading e−γt Explosive blast loading (29) where t1 = 0.006 s, γ = 330 s−1 and q0 = 3.448 MPa. Table 4. Properties of material Material E1(GPa) E2 (GPa) G12 (GPa) G13 (GPa) G23 (GPa) ν12 ρ (kg/m3) I 172.369 6.895 3.448 3.448 3.448 0.25 1603.03 II 172.369 6.895 3.448 3.448 1.379 0.25 1603.03 III 131.69 8.55 6.67 6.67 6.67 0.3 1610 Fig. 3 shows the time histories of central deflection of the plate under various dy- namic loadings. The obtained results of present solution using IGA are compared with those obtained by Wang et al. [7] using the strip element method (SEM). As expected, the effectiveness of this work is fully believable when profiles relatively coincide with Wang et al.’s solutions. 276 Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan 9 a) step loading b) Triangular loading c) sine loading d) explosive blast loading Fig. 3. Variation of the center deflection as a function of time for a (00/900/00) square laminated composite plate subjected to various dynamic loadings 4.2.3 A circular four-layer plate [45 0 /-45 0 /-45 0 /45 0 ] Finally, to increase lively for numerical examples and obtain the desired effect, we consider a [450/-450/-450/450] circular plate with fully clamped (CCCC) boundary condition as shown Fig. 6a. Material parameter III is used. The plate is also subjected to a conventional blast load as given in Section 4.2.2. The circular plate has the radius to thickness ratio is 10 (R/h = 10). A rational quadratic basis is enough to model exactly the circular geometry. Coarse mesh and control net of the plate with respect to quadratic and cubic elements are illustrated in Fig. 7. Time step for transient analysis is chosen  t = 0.4ms. The plate is meshed with 13x13 NURBS cubic elements as shown Fig. 6b. Fig. 8 illustrates the profile of displacements versus time at the center of the circular plate subjected to conventional blast load. Obtained results are compared with solutions from ANSYS 13 which using SHELL 181 elements. It can be seen that the present solutions are in good agreement with the solutions from ANSYS software. (a) Step loading a) step loading b) Triangular loading c) sine loading d) explosive blast loading Fig. 3. Variation of the center deflection as a function of time for a (00/900/00) square laminated composite plate subjected to various dynamic loadings 4.2.3 A circular four-layer plate [45 0 /-45 0 /-45 0 /45 0 ] Finally, to increase lively for numerical examples and obtain the desired effect, we consider a [450/-450/-450/450] circular plate with fully clamped (CCCC) boundary condition as shown Fig. 6a. Material parameter III is used. The plate is also subjected to a conventional blast load as given in Section 4.2.2. The circular plate has the radius to thickness ratio is 10 (R/h = 10). A rational quadratic basis is enough to model exactly the circular geometry. Coarse mesh and control net of the plate with respect to quadratic and cubic elements are illustrated in Fig. 7. Time step for transient analysis is chosen  t = 0.4ms. The plate is meshed with 13x13 NURBS cubic elements as shown Fig. 6b. Fig. 8 illustrates the profile of displacements versus time at the center of the circular plate subjected to conventional blast load. Obtained results are compared with solutions from ANSYS 13 which using SHELL 181 elements. It can be seen that the present solutions are in good agreement with the solutions from ANSYS software. (b) Triangular loading 9 a) step loading b) Triangular loading c) sine loading d) explosive blast loading Fig. 3. Variation of the center deflection as a function of time for a (00/900/00) square laminated composite plate subjected to various dynamic loadings 4.2.3 A circular four-layer plate [45 0 /-45 0 /-45 0 /45 0 ] Finally, to increase lively for numerical examples and obtain the desired effect, we consider a [450/-450/-450/450] circular plate with fully clamped (CCCC) boundary condition as shown Fig. 6a. Material parameter III is used. The plate is also subjected to a conventional blast load as given in Section 4.2.2. The circular plate has the radius to thickness ratio is 10 (R/h = 10). A rational quadratic basis is enough to model exactly the circular geometry. Coarse mesh and control net of the plate with respect to quadratic and cubic elements are illustrated in Fig. 7. Time step for transient analysis is hos  t = 0.4ms. The plate is m hed with 13x13 NURBS c bic elements as shown Fig. 6b. Fi . 8 illustrates the profile of displacem nts versus time at the center of the circular plat subjected to c nventional blast load. Obt ined results are compared with solutions from ANSYS 13 which using SHELL 181 elements. It can be seen that the present solutions are in good agreement with the solutions from ANSYS software. (c) S ne l ading 9 a) step loading b) Triangular loading c) sine loading d) explosive blast loading Fig. 3. Variation of the center deflection as a function of time for a (00/900/00) square laminated composite plate subjected to various dynamic loadings 4.2.3 circular four-layer plate [45 0 /-45 0 /-45 0 /45 0 ] Final y, to increase lively for numerical examples and obtain the desired effect, we consider a [450/-450/-450/450] circular plate with fully clamped (C ) boundary condit on as shown Fig. 6a. aterial parameter I is used. The plate is also subjected to a conventional blast load as given in Section 4.2.2. The circular plate has the radius to thicknes ratio is 10 (R/h = 10). A rational quadratic basis is enough to model exactly the circular geometry. Coarse mesh and control net of the plate with respect to quadratic and cubic elements are illustrated in Fig. 7. Time step for transient analysi is chosen  t = 0.4ms. The pla e is meshed with 13x13 NURBS cubic elements as hown Fig. 6b. Fig. 8 il ustrates the profile of displacements versus time at he center of the circular plate subjected to conventional blast load. Obtain d results are compared with solutions from ANSYS 13 which using S ELL 181 elements. It can be se n that the present solutions are in go d agreement with the solutions from ANSYS software. (d) Explosive blast loading Fig. 3. Variation of the c nter deflec ion as a fun ion of time for (00/900/00) square laminated composite plate subjected to various dynamic loadings Second, a fully simply s pported thre -l yer square laminat d plate s rte as 00/900 00] is also onsider d. Material II is used. The length and thickne of the plates are assumed to be a = 5h and h = 0.1524 , respectively. As above example, the plat is al o subjecte to si usoid lly distributed transv rse load (with q0 68.9476 MPa). The displacement at the center of plat is also studied. Khdeir and Reddy [21] orig nally inves- tiga ed this benchmark olution. Fig. 4 hows variation of the displacement at th c nter of plate as a function under various dynamic loadings. The present solutions based on IGA and TSDT ar compared with exact solution of Khdeir and Reddy [21] using HSDT. As observed in Fig. 4, the profiles are relatively accurate, the error estimate is very small and approvable when comparing with exact solution. Transient analysis of laminated composite plates using NURBS-based isogeometric analysis 277 10 a) step loading b) Triangular loading c) sine loading d) explosive blast loading Fig. 4. Central deflection versus time for a [00/900/00] square laminated plate subjected to various dynamic loadings Fig. 5. The time history of the center deflection of the [300/-300/-300/300] fully clamped laminated plate. (a) Step loading 10 i g b) Triangular loading c) sine loading d) explosive blast loading Fig. 4. entral deflection versus time for a [00/900/00] square laminated plate subjected to various dynamic loadings Fig. 5. The time history of the center deflection of the [300/-300/-300/300] fully clamped laminated plate. (b) Triangular loading 10 a) step loading b) Triangular loading c) sin oading d) explosive blast loading Fig. 4. Central deflection versus time for a [00/900/00] square laminated plate subjected to various dynamic loadings Fig. 5. The time history of the center deflection of the [300/-300/-300/300] fully clamped laminated plate. (c) Sin loading 10 a) step loading b) Triangular loading c) sine loading d) explosive blast loading Fig. 4. Central deflection versus time for a [00/900/00] square laminated plate subjected to various dynamic loadings Fig. 5. The time history of the cent r deflection f the [300/-300/-300/300] fully clamped laminated plate. (d) Explosive blast loading Fig. 4. Central deflection versus time for a [00/900/00] square laminated plate subjected to various dynamic loadings 4.2.2. A four-layer square plate [300/− 300/− 300/300] A fully clamped four-layer angle-ply square laminated plate with symmetrically stacking sequences [300/− 300/− 300/300] is considered. Material III is used. The length to thickness ratio of the plate is assumed to be a/h = 50. The plate is also subjected to a transverse load which is uniformly distributed over the plate and is called conventional blast loading [7]. q(x, y, t) = q0(1− t t2 )e−α1t/t2 , (30) in which q0 = 68.9476 KPa, t2 = 0.004 s, α1 = 1.98. 278 Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan 10 a) step loading b) Triangular loading c) sine loading d) explosive blast loading Fig. 4. Central deflection versus time for a [00/900/00] square laminated plate subjected to various dynamic loadings Fig. 5. The time history of the center deflection of the [300/-300/-300/300] fully clamped laminated plate. Fig. 5. The time history of the center deflection of the [300/− 300/− 300/300] fully clamped laminated plate The time history of the deflection at the center of the four-layer fully clamped (CCCC) laminated plate is investigated, as shown in Fig. 5. The results are compared with the solutions of Wang et al. [7]. From Fig. 5, the present results match well with the reference solutions. 4.2.3. A circular four-layer plate [450/− 450/− 450/450] Finally, to increase lively for numerical examples and obtain the desired effect, we consider a [450/ − 450/ − 450/450] circular plate with fully clamped (CCCC) boundary condition as shown Fig. 6a. Material parameter III is used. The plate is also subjected to a conventional blast load as given in Section 4.2.2. The circular plate has the radius to thickness ratio is 10 (R/h = 10). A rational quadratic basis is enough to model exactly the circular geometry. Coarse mesh and control net of the plate with respect to quadratic and cubic elements are illustrated in Fig. 7. Time step for transient analysis is chosen ∆t = 0.4 ms. The plate is meshed with 13× 13 NURBS cubic elements as shown Fig. 6b. Fig. 8 illustrates the profile of displacements versus time at the center of the circular plate 11 Fig. 6. The circular plate: (a) geometry and (b) mesh based on 13x13 cubic elements. Fig. 7. Coarse mesh and control points of a circular plate with various degrees: a) p=2 and b) p=3. 5. CONCLUSIONS Isogeometric analysis combined with TSDT to analyze the transient of laminated composite plates is first studied. The displacement field is generally defined and is derived from CPT. The Newmark time-integration algorithm was chosen to approximate the ordinary differential equations in time. We have successfully extended an application of the NURBS-based isogeometric finite element approach to transient analysis for laminated composite plates as this work. IGA is the effectively numerical method. It has expressed well its role in solving the problems with just few elements especially curved geometry as circle. The calculation of these problems has been done very fast. It not only helps to save costs but also increases the accuracy of solutions. The numerical results agreed well with those of available references and exact solution, and hence illustrated the accuracy and effectiveness of the present method. Fig. 8. The deflection at the center of the [450/-450/-450/450] circular laminated plate subjected to a conventional blast load. (a) 11 Fig. 6. The circular plate: (a) geometry and (b) mesh based on 13x13 cubic elements. Fig. 7. Coarse mesh and control points of a circular plate with various degrees: a) p=2 and b) p=3. 5. CONCLUSIONS Isogeometric analysis combined with TSDT to analyze the transient of la inated c mposite plates is first studied. The displacement field is generally defined and is derived from CPT. The Newmark time-integration algorithm was chosen to approximate the ordinary differential equations in time. We have successfully extended an application of the NURBS-based isogeometric finite element approach to transient analysis for laminated composite plates as this work. IGA is the effectively numerical method. It has expressed well its role in solving the problems with just few elements especially curved geometry as circle. The calculation of these problems has been done very fast. It not only helps to save costs but also increases the accuracy of solutions. The numerical results agreed well with those of available references and exact solution, and hence illustrated the accuracy and effectiveness of the present method. Fig. 8. The deflection at the center of the [450/-450/-450/450] circular laminated plate subjected to a conventional blast load. (b) Fig. 6. The circular plate: (a) geometry and (b) mesh based on 13× 13 cubic elements Transient analysis of laminated composite plates using NURBS-based isogeometric analysis 279 subjected to conventional blast load. Obtained results are compared with solutions from ANSYS 13 which using SHELL 181 elements. It can be seen that the present solutions are in good agreement with the solutions from ANSYS software. 11 Fig. 6. The circular plate: (a) geometry and (b) mesh based on 13x13 cubic elements. Fig. 7. Coarse mesh and control points of a circular plate with various degrees: a) p=2 and b) p=3. 5. CONCLUSIONS Isogeometric analysis combined with TSDT to analyze the transient of laminated composite plates is first studied. The displacement field is generally defined and is derived from CPT. The Newmark time-integration algorithm was chosen to approximate the ordinary differential equations in time. We have successfully extended an application of the NURBS-based isogeometric finite element approach to transient analysis for laminated composite plates as this work. IGA is the effectively numerical method. It has expressed well its role in solving the problems with just few elements especially curved geometry as circle. The calculation of these problems has been done very fast. It not only helps to save costs but also increases the accuracy of solutions. The numerical results agreed well with those of available references and exact solution, and hence illustrated the accuracy and effectiveness of the present method. Fig. 8. The deflection at the center of the [450/-450/-450/450] circular laminated plate subjected to a conventional blast load. (a) p = 2 11 Fig. 6. The circular plate: (a) geometry and (b) mesh based on 13x13 cubic elements. Fig. 7. Coarse mesh nd control points of a ci cular plate with various degrees: a) p=2 and b) p=3. 5. CONCLUSIONS Isogeometric analysis combined with TSDT to analyze the transien of laminated composite plates is first stud ed. The displacement field is generally defined and is derived from CPT. The Newmark time-int gration algorithm was chosen to approximate the ordinary differential equations in time. We have su cessfully extended an pplication of he NURBS-based isogeometric finite element approach to transient analysis for laminated composite plates as his work. IGA is the effectively numerical meth d. It has expressed well it role in olving the problems with just few elements especially curved geometry as circle. The calculatio of thes problems has been done very fast. It not only helps to save costs but also increases the accu acy of solutions. The numerical r sults agreed well with those of available references and exac solution, and hence illustrat d the accuracy and effectiveness of the present method. Fig. 8. The deflection at the cente of the [450/-450/-450/450] circular laminated plate subjected to a conventional blast load. (b) p = 3 Fig. 7. Coarse mesh and control points of a circular plate with various degrees 11 Fig. 6. The circular plate: (a) geometry and (b) mesh based on 13x13 cubic elements. Fig. 7. Coarse mesh and control points of a circular plate with various degrees: a) p=2 and b) p=3. 5. CONCLUSIONS Isogeometric analysis combined with TSDT to analyze the transient of laminated composite plates is first studied. The displacement field is generally defined and is derived from CPT. The Newmark time-integration algorithm was chosen to approximate the ordinary differential equations in time. We have successfully extended an application of the NURBS-based isogeometric finite element approach to transient analysis for laminated composite plates as this work. IGA is the effectively numerical method. It has expressed well its role in solving the problems with just few elements especially curved geometry as circle. The calculation of these problems has been done very fast. It not only helps to save costs but also increases the accuracy of solutions. The numerical results agreed well with those of available references and exact solution, and hence illustrated the accuracy and effectiveness of the present method. Fig. 8. The deflection at the center of the [450/-450/-450/450] circular laminated plate subjected to a conventional blast load. Fig. 8. The deflection at the center of the [450/− 450/− 450/450] circular laminated plate subjected to a conventional blast load 5. CONCLUSIONS Isogeometric analysis combined with TSDT to analyze the transient of laminated composite plates is first studied. The displacement field is generally defined and is de- rived from CPT. The Newmark time-integration algorithm was chosen to approximate the ordinary differential equations in time. We have successfully extended an application of the NURBS-based isogeometric finite element approach to analyze dynamic response for laminated composite plates as this work. IGA is the effectively numerical method. It has expressed well its role in solving the problems with just few elements especially curved 280 Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan geometry as circle. The calculation of these problems has been done very fast. It not only helps to save costs but also increases the accuracy of solutions. The numerical results agreed well with those of available references and exact solution, and hence illustrated the accuracy and effectiveness of the present method. ACKNOWLEDGMENT This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.02-2012.17. The support is gratefully acknowledged. REFERENCES [1] H. Reismann. Forced motion of elastic plates. Journal of Applied Mechanics, 35, (3), (1968), pp. 510–515. [2] H. Reismann and Y. C. Lee. Forced motion of rectangular plates. Developments in Theoretical and Applied Mechanics, 4, (1968), pp. 3–18. [3] T. Rock and E. Hinton. Free vibration and transient response of thick and thin plates using the finite element method. Earthquake Engineering & Structural Dynamics, 3, (1), (1974), pp. 51–63. [4] H. U. Akay. Dynamic large deflection analysis of plates using mixed finite elements. Computers & Structures, 11, (1), (1980), pp. 1–11. [5] J. N. Reddy. Dynamic (transient) analysis of layered anisotropic composite-material plates. International Journal for Numerical Methods in Engineering, 19, (2), (1983), pp. 237–255. [6] Mallikarjuna and T. Kant. Dynamics of laminated composite plates with a higher order theory and finite element discretization. Journal of Sound and Vibration, 126, (3), (1988), pp. 463– 475. [7] Y. Wang, K. Lam, and G. Liu. A strip element method for the transient analysis of symmetric laminated plates. International Journal of Solids and Structures, 38, (2), (2001), pp. 241–259. [8] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. Isogeometric analysis: Cad, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194, (39), (2005), pp. 4135–4195. [9] C. H. Thai, A. J. M. Ferreira, S. P. A. Bordas, T. Rabczuk, and H. Nguyen-Xuan. Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. European Journal of Mechanics-A/Solids, 43, (2014), pp. 89–108. [10] N. Valizadeh, S. S. Ghorashi, H. Yousefi, T. Q. Bui, and T. Rabczuk. Transient analysis of laminated composite plates using isogeometric analysis. In Proceedings of the Eighth Interna- tional Conference on Engineering Computational Technology, Civil-Comp Press, Stirlingshire, Scotland, (2012). [11] J. N. Reddy. A simple higher-order theory for laminated composite plates. Journal of Applied Mechanics, 51, (4), (1984), pp. 745–752. [12] J. N. Reddy. Mechanics of laminated composite plates and shells: Theory and analysis. CRC press, (2004). [13] N. M. Newmark. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 85, (3), (1959), pp. 67–94. [14] A. A. Khdeir and L. Librescu. Analysis of symmetric cross-ply laminated elastic plates using a higher-order theory: Part II-buckling and free vibration. Composite Structures, 9, (4), (1988), pp. 259–277. Transient analysis of laminated composite plates using NURBS-based isogeometric analysis 281 [15] J. N. Reddy. Mechanics of laminated composite plates: Theory and analysis. CRC press Boca Raton, (1997). [16] K. M. Liew, Y. Q. Huang, and J. N. Reddy. Vibration analysis of symmetrically laminated plates based on FSDT using the moving least squares differential quadrature method. Com- puter Methods in Applied Mechanics and Engineering, 192, (19), (2003), pp. 2203–2222. [17] A. J. M. Ferreira. A formulation of the multiquadric radial basis function method for the analysis of laminated composite plates. Composite Structures, 59, (3), (2003), pp. 385–392. [18] A. J. M. Ferreira, L. Castro, and S. Bertoluzza. A high order collocation method for the static and vibration analysis of composite plates using a first-order theory. Composite Structures, 89, (3), (2009), pp. 424–432. [19] W. Zhen and C. Wanji. Free vibration of laminated composite and sandwich plates using global-local higher-order theory. Journal of Sound and Vibration, 298, (1), (2006), pp. 333– 349. [20] H. Matsunaga. Vibration and stability of cross-ply laminated composite plates according to a global higher-order plate theory. Composite Structures, 48, (4), (2000), pp. 231–244. [21] A. A. Khdeir and J. N. Reddy. Exact solutions for the transient response of symmetric cross- ply laminates using a higher-order plate theory. Composites Science and Technology, 34, (3), (1989), pp. 205–224. VIETNAM ACADEMY OF SCIENCE AND TECHNOLOGY VIETNAM JOURNAL OF MECHANICS VOLUME 36, N. 4, 2014 CONTENTS Pages 1. N. T. Khiem, P. T. Hang, Spectral analysis of multiple cracked beam subjected to moving load. 245 2. Dao Van Dung, Vu Hoai Nam, An analytical approach to analyze nonlin- ear dynamic response of eccentrically stiffened functionally graded circular cylindrical shells subjected to time dependent axial compression and external pressure. Part 2: Numerical results and discussion. 255 3. Lieu B. Nguyen, Chien H. Thai, Ngon T. Dang, H. Nguyen-Xuan, Tran- sient analysis of laminated composite plates using NURBS-based isogeometric analysis. 267 4. Tran Xuan Bo, Pham Tat Thang, Do Thanh Cong, Ngo Sy Loc, Experimental investigation of friction behavior in pre-sliding regime for pneumatic cylinder 283 5. Vu Thi Thuy Anh, Dao Huy Bich, Nguyen Dinh Duc, Nonlinear post-buckling of thin FGM annular spherical shells under mechanical loads and resting on elastic foundations. 291 6. N. D. Anh, N. N. Linh, A weighted dual criterion for stochastic equivalent linearization method using piecewise linear functions. 307