Dynamic instability of curved variable angle tow composite panel under axial compression

Dynamic instability of curved variable angle tow composite panel under axial compression

Thin-Walled Structures 138 (2019) 302–312 Contents lists available at ScienceDirect Thin-Walled Structures journal homepage: www.elsevier.com/locate...

3MB Sizes 0 Downloads 44 Views

Thin-Walled Structures 138 (2019) 302–312

Contents lists available at ScienceDirect

Thin-Walled Structures journal homepage: www.elsevier.com/locate/tws

Full length article

Dynamic instability of curved variable angle tow composite panel under axial compression

T



Surya Samukhama, Gangadharan Rajua, , C.P. Vyasarayania, Paul M. Weaverb a b

Department of Mechanical and Aerospace Engineering, Indian Institute of Technology Hyderabad, Telangana 502205, India Bernal Institute, School of Engineering, University of Limerick, Castletroy, Ireland

A R T I C LE I N FO

A B S T R A C T

Keywords: Dynamic instability Variable angle tow composite Bolotin's approach Generalized differential integral quadrature method Curved panel

Variable angle tow (VAT) composites have demonstrated better performance in buckling and post-buckling over straight fiber composites based on the mechanics of load redistribution from critical regions to supported edges. In this work, the dynamic instability behavior of a curved VAT composite panel subjected to periodic axial compression load is investigated. The governing energy functional of a curved symmetric VAT panel under external loading is derived using Donnell's shallow shell theory. Later, the discretized equations of motion are derived using the Rayleigh–Ritz method combined with the generalized differential integral quadrature method (GDIQM). Initially, the pre-buckling problem is solved by applying a uniform compression load to compute the stress resultant distribution which is used to evaluate the buckling load of the curved VAT panel. Subsequently, the dynamic/parametric instability region of a curved VAT panel subjected to periodic axial compression load is determined using Bolotin's first-order approximation. Then, the dynamic instability performance is evaluated for a curved VAT panel with linear fiber angle distribution and compared with straight fiber laminates. Finally, the influence of fiber angle orientation, the radius of curvature, aspect ratio and plate boundary conditions on the dynamic instability of VAT panel is presented.

1. Introduction The emergence of variable angle tow composite fabrication technology in recent years has provided more options to tailor structural stiffness distributions that enhance mechanical behavior under different loading conditions. The mechanics underpinning the improvement of buckling and post-buckling performance of a VAT laminate is primarily attributed to load redistribution from critical regions to supported edges of the laminate. In the works of Hyer and Lee [1] and Hyer and Charette [1], showed that the VAT panels can improve the stress concentration around holes by arranging the fibers in the direction of critical load paths. Gurdal and Olmedo [2] introduced a linear fiber angle variation along the length and studied the in-plane response of composite laminates with variable stiffness. Later, Gurdal et al. [3] studied the effect of stiffness variation on buckling response of VAT composite laminates using Rayleigh–Ritz method. In their studies, they have demonstrated the load re-distribution phenomenon is responsible for the improvement of buckling performance of the VAT composites. After that numerous works in the literature have reported the benefits of VAT composites over straight fiber composites [4–9]. In the recent times, the advancement in the automation of composite manufacturing enables



the fabrication of VAT panels with greater reliability and minimum defects [10,11]. In aircraft structures, the primary load carrying parts such as fuselage and wing are subjected to dynamic loads under different flight conditions. These structures when subjected to time-periodic in-plane compression loads, undergo parametric resonance for certain combinations of excitation frequency and amplitude of axial load. As a result, instability may occur below the critical load of the structure over a range of excitation frequencies. Therefore, understanding the dynamic instability characteristics of a structure are of great practical important from the theoretical and design perspective. A comprehensive study on the dynamic/parametric instability behavior of various elastic systems has been carried out by Bolotin [12]. Many researchers have used various modeling approaches to study the dynamic instability of structures using Bolotin's method. Srinivasan and Chellapandi [13] used the finite strip method (FSM) based on classical laminated plate theory to perform the dynamic instability analysis of laminated plates under periodic in-plane load. Moorthy et al. [14], Chattopadhyay and Radu [15] and Dey and Singha [16] carried out a similar investigation using the finite element method (FEM) to approximate the instability regions for moderately thick composite plates. Sahu and Datta [17]

Corresponding author. E-mail address: [email protected] (G. Raju).

https://doi.org/10.1016/j.tws.2019.02.015 Received 26 October 2018; Received in revised form 7 January 2019; Accepted 13 February 2019 0263-8231/ © 2019 Elsevier Ltd. All rights reserved.

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

used FEM combined with Bolotin's approach for the parametric instability of composite laminates subjected to non-uniform periodic loads. Cederbaum [18] used method of multiple scales to determine the dynamic instability of composite panels subjected to in-plane periodic loads. Wang and Dawe [19] implemented a B-spline finite strip method based on first order shear deformation theory for the dynamic instability analysis of rectangular and prismatic composite plate structures. Chen et al. [20] and Samukham et al. [21] investigated the dynamic instability of a VAT composite laminate with and without delamination, respectively. In both works, considerable improvement in dynamic stability of VAT laminate over straight fiber laminates was demonstrated. Recently, Khalafi and Fazilati [22] studied the effect of geometry and fiber angle orientation on parametric instability of VAT quadrilateral plates using iso-geometric analysis with finite element formulation. All of the previous analyses were performed on flat plates, whereas, in modern aerospace, civil, mechanical and naval engineering structures curved panels are extensively used for load bearing applications. Curved panels can exhibit greater structural efficiency because of their ability to sustain high compressive loads. The improvement in their buckling performance is mainly attributed to the curvature of the panel and to their ability to develop membrane loads. Ganapathi et al. [23] studied the effect of curvature and aspect ratio on dynamic instability of a uniformly loaded thick composite cylindrical panels using the finite element method. They demonstrated the effectiveness of a nine-node shear flexible shell element in computing the dynamic instability region of a curved composite panel. Sahu and Datta [24] used finite element analysis to investigate the dynamic stability of single and doubly curved composite panels including spherical, elliptic, and hyperbolic paraboloids. The dynamic stability of thin composite cylindrical shells under combined static and periodic axial force was investigated by Ng et al. [25] and Lam and Ng [26] using Love's theory of thin shells. Many works on dynamic instability of composite laminates are based on finite element models and suffer from the issue of shear locking for thin laminates. To overcome this problem, we have used the recently developed generalized differential integral quadrature method to solve the dynamic instability of VAT laminate. From the literature, it is evident that a large amount of work has been conducted on the parametric instability of straight fiber curved panels. However, it appears no work has been reported on the parametric instability of curved VAT panels. In this work, dynamic instability of a symmetric curved VAT panel subjected to a periodic uniaxial compression is studied. The generalized differential-integral quadrature method (GDIQM) is used to solve the integral energy expressions and Bolotin's first-order approximation is used to determine the dynamic instability boundary of a curved VAT panel. The effect of tow steering, radius of curvature, aspect ratio and plate boundary conditions on buckling and dynamic instability behavior of curved VAT panel under periodic axial compression load is investigated, and compared with straight fiber composite panels. The paper is organized as follows: In Section 2, we discuss the modeling of a curved panel with VAT laminate. In Section 3, the details of GDIQM is explained. In Section 4 Bolotin's first-order approximation for obtaining stability boundary of a VAT panel under periodic in-plane compression load is discussed. We validate our model in Section 5 with results from a commercial finite element software (ABAQUS) and results from the literature. Furthermore, the influence of fiber angle orientation, radius of curvature and aspect ratio on the dynamic instability of a curved VAT panel is presented. Finally, we conclude our observations in Section 6.

Fig. 1. Schematic representation of a curved panel with VAT laminate.

linear fiber angle variation over the lamina, Lagrange polynomials [2] can be used as:

ψ (x ) = ϕ +

2(T1 − T0) |x| + T0, a

(1)

where, T0 is the fiber angle at the center of the laminate (x = 0) and T1 is the fiber angle at the laminate ends (x = ± a/2) and ϕ is the angle of rotation of the fiber path. We assume the thickness of the laminate to be small, therefore, the shear deformation of the shell's section in x − z and y − z planes is neglected. Using Donnell's shallow shell kinematic equations [27] and neglecting higher order strain terms, the mid-plane strains (ϵ ) and curvatures (κ ) of the laminate can be written as: T

∂u ∂v ∂v ⎞ ⎤ w ∂u + ⎞⎛ + ϵ =⎡ ⎛ , ⎢ ∂x ∂y ∂ ∂ R y x ⎠⎥ ⎝ ⎠ ⎝ ⎦ ⎣ ⎜

⎟⎜



T

∂ 2w ∂ 2w ∂ 2w ⎤ κ = −⎡ 2 22 , ⎢ ⎣ ∂x ∂y ∂x ∂y ⎥ ⎦

(2)

where, u, v, w are the mid-plane displacement fields corresponding to x , y and z directions, respectively. The variation of fiber angle orientation over the laminate results in the following constitutive relationship:

¯ (x , y ) B ¯ (x , y ) ⎤ ⎧ ϵ ⎫ ¯ ⎫ ⎡A ⎧N = , ¯ ⎬ ⎢B ¯ T (x , y ) D ¯ (x , y ) ⎥ ⎨ κ ⎬ ⎨ M ⎩ ⎭ ⎣ ⎦⎩ ⎭

(3)

¯ = [N¯x N¯ y N¯xy]T is the membrane stress resultant vector and where, N ¯ = [M ¯x M ¯yM ¯ xy]T is the bending moment resultant vector. Since the M fiber angle orientation varies with respect to position over the lamina, ¯ , coupling stiffness matrix B ¯ and the the in-plane stiffness matrix A ¯ are all functions of spatial coordinates (x , y ) . bending stiffness matrix D The strain energy [21] of a VAT panel subjected to uni-axial compression is given by: U=

1 2

∫ [ϵT A¯ ϵ + κT D¯ κ] dA + ∫ [ϵLT T¯ ϵL] d ,

(4)

¯ and ϵL are defined as follows: where, T N¯ N¯xy ⎤ ¯ =⎡ x T ⎢ N¯xy N¯ y ⎥ ⎣ ⎦

(5) T

∂w ∂w ⎤ ϵL = ⎡ . ⎢ ⎣ ∂x ∂y ⎥ ⎦

(6)

In this work, we combine the Rayleigh–Ritz method with GDIQM and apply them to the strain energy functional of the VAT panel Eq. (4) and solve for the dynamic instability. Unlike straight fiber composites, the varying stiffness properties of the VAT panel induces a non-uniform stress resultant distribution across the plane of the panel under pure compression. Therefore, the pre-buckling analysis of a curved VAT

2. Modeling of VAT laminate A curved panel with VAT laminate considered is shown in Fig. 1. In VAT laminates, the fiber angle orientation continuously changes with respect to x − y coordinates in the plane of the laminate. To describe a 303

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

panel is carried out initially by applying uniform compression load to obtain the stress resultant distribution. Next, the obtained stress resultant distribution is used to calculate the geometric stiffness matrix K G . From the obtained K G matrix the dynamic instability region of a VAT panel is determined from the following equations of motion:

matrices along the x and y directions, respectively. In the same way, higher-order derivative weighting coefficient matrices can be derived on the basis of Eq. (8) (i.e., by matrix multiplication) or by recursive formulas [32]. For example, second-order partial derivatives can be written as:

Mu ¨ + Ku − (α 0 Pcr + α1 Pcr cos(θt )) K G u = 0,

∂2 f = T xxf, ∂x 2

(7)

where, K represents the structural stiffness matrix, M is the mass matrix, Pcr is the critical buckling load of the VAT panel, α 0 and α1 are static and dynamic load parameters, respectively and θ is the excitation frequency. In the next sections, we discuss the application of GDIQM and the Rayleigh–Ritz method for obtaining the explicit expressions for M, K , and K G in the equation of motion (Eq. (7)).

T xx ,

π (j − 1) ⎞ ⎞ ⎛ b ⎜1 − cos ⎜⎛ ⎟ ⎟ ⎝ (n y − 1) ⎠ ⎠ ⎝ yi = . 2

T ikx fk with k = 1, …, n x n y .

∫ g (x ) dx ≈ P xg + c

∂ f = T yf , ∂y

with P x = (T x)+,

(12)

()+

is the Moore-Penrose pseudo-inverse operator, the vector g where contains the value of the function at grid points and c is the vector of integration constants. Here P x has the same dimension as T x . In order to approximate one dimensional definite integrals, we define the following vectors j x and j y as: y y x y j x = P(xnx , i) − P(1, i) , i = 1, …, n x and j = P(ny, k ) − P(1, k ) , k = 1, …, n y ,

(8)

(13)

Here, f (x , y ) is discretized into fi for the grid (x , y )i and n x , n y represents the number of grid points along x and y directions, respectively (see Fig. 2). The matrix T x is the weighting coefficient matrix of firstorder derivative. By converting the grid point variable fi into a vector form f , first-order partial derivative can be written in matrix form as:

∂ f = T xf, ∂x

(11)

The schematic representation of a discretized Chebyshev grid of a curved panel along with the GDIQM numbering system is shown in Fig. 2. In GDIQM, the integral operations in the strain energy functional Eq. (4) are replaced by an integral-quadrature matrix operator [32]. In this work, the pseudo-inverse operator is used [35] to compute the integral weighting coefficient matrices. For a one-dimensional function g (x ) , the weighting coefficient matrix for an arbitrary grid spacing can be written as:

nx × ny i=1

(10)

T xy

T yy

π (i − 1) ⎞ ⎞ ⎟⎟ a ⎛⎜1 − cos ⎜⎛ ⎝ (n x − 1) ⎠ ⎠ xi = ⎝ , 2

In the literature, buckling and post-buckling problems of VAT laminates have been successfully solved with the differential quadrature method [4,28] (DQM). The DQM approach is based on the approximation of a derivative using a weighted linear sum of the function values at discretized grid points in the domain. It should be noted that on application of DQM discretization, the partial differential equation (PDE) governing the plate behavior is converted into a system of algebraic equations. The numerical incorporation of strong form of plate boundary conditions along the boundary grid points in DQM is non trivial and challenging. To overcome this difficulty, the problem is expressed in weak form using the strain energy functional, and then the Rayleigh–Ritz combined with GDIQM is used to get the equilibrium equations. However, one of the major limitation of GDIQM is that it cannot be applied to the structures with geometric discontinuities. White et al. [29] introduced the GDIQM approach to solve the integrodifferential post-buckling equations of a VAT laminate. Raju et al. [30] applied a perturbation based asymptotic numerical method coupled with GDIQM to perform the post-buckling of VAT laminates. In this work, GDIQM is used to solve the problems of pre-buckling, buckling and dynamic instability of a curved VAT panel subjected to a periodic axial compression load. In DQM [31], the partial derivative of a function f (x , y ) can be written as:



∂2 f = T xyf, ∂x ∂y

and are the second-order derivative weighting where, coefficient matrices. Numerical values of the weighting coefficients Tijx are dependent on DQ basis function and spacing of the grid points. For plate bending and buckling problems, the DQM solutions converged quickly when Lagrange polynomials were used as basis functions along with Chebyshev grid spacing [4,33,34]. The Chebyshev grid coordinates are:

3. Generalized differential-integral quadrature method (GDIQM)

∂ f = ∂x i

∂2 f = T yyf, ∂y 2

jx

jy

where, the vectors and have dimensions 1 × n x and 1 × n y , respectively. Using the vectors j x and j y the definite integrals of a function g (x , y ) along the lines of constant y or x can be constructed as: T a ∫0 ⎡⎢g ⎜⎛x , y1⎟⎞, g ⎜⎛x , y2 ⎟⎞, …, g ⎜⎛x , yny ⎟⎞ ⎤⎥ dx ≈ J xg

⎠ ⎝ ⎠ ⎝ ⎠⎦ ⎣ ⎝ T dx ≈ J yg , [ g ( x , y ), g ( x , y ), … , g ( x , y )] ∫0 nx 1 2

(9)

a

where, T x and T y are the first-order derivative weighting coefficient

Jx

(14)

Jy

where, the matrices and contains the elements of vectors j x and j y arranged according to the grid-point numbering scheme. When the system defined in Fig. 2 is used, then the matrices can be computed by the matrix Kronecker product [36], denoted by ⊗. Therefore,

J x = j x ⊗ I y and J y = Ix ⊗ j y ,

Ix

nx

Iy

 ny

(15)

Jx

ny × nx ny

∈ ∈ ∈ and are identity matrices and , where J y ∈ nx × nx ny are two dimensional integral weighting matrices. For evaluating, two-dimensional integrals, Eq. (14) is written as: a

∫0 ∫0 Fig. 2. A two-dimensional Chebyschev grid plotted on the surface of a curved panel with the grid-point numbering system.

b

g ⎛⎜x , y ⎞⎟ dxdy ≈ ⎛⎜j y J x⎞⎟ g = s xyg . ⎝ ⎠ ⎝ ⎠

(16)

The integral of a product of two functions f (x , y ) and g (x , y ) can be calculated approximately by constructing a matrix S = diag(s xy) : 304

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

a

∫0 ∫0

b

f ⎛⎜x , y ⎞⎟ g ⎛⎜x , y ⎞⎟ dxdy ≈ f T Sg . ⎝ ⎠ ⎝ ⎠

(17)

In the next section, the GDIQM approach is applied to Eq. (4) using the weighting coefficient matrices T x , T y , T xx , T yy , T xy and S . 3.0.1. Pre-buckling analysis On applying Rayleigh–Ritz to the first term in the strain energy functional (Eq. (4)), the system of equations that governs the prebuckling behavior of a VAT panel can be written as:

⎡ Kuu Kuv Kuw ⎤ ⎧ u ⎫ ⎧ Fu ⎫ ⎢ Kvu Kvv Kvw ⎥ v = Fv , w⎬ ⎨ ⎬ ⎢ Kwu Kwv Kww ⎥ ⎨ ⎦ ⎩ ⎭ ⎩ Fw ⎭ ⎣

(18) Fig. 3. Schematic representation of in-plane boundary conditions for prebuckling analysis of a curved panel with VAT laminate.

and can be written in compact form as: (19)

Ks u = F,

where, u represents the pre-buckling solution and F represents the forces applied along the x, y and z directions, respectively. The terms of Ks matrix in Eq. (18) are expressed as follows:

using the equations uslave − ucontrol = 0 , which is similar to a multi-point constraint (MPC) approach used in the commercial finite element analysis. These constraint equations can be written in matrix form as K c u = 0 . The constraint matrix, K c is added to the structural stiffness matrix Ks to obtain the final stiffness matrix K .

∂′ ∂ ∂′ ∂ ∂′ ∂ ∂′ ∂ + + + Kuu = ∫ ⎡ A11 A16 A16 A66 ⎤ d  ⎢ ∂x ∂y ∂x ∂x ∂y ∂y ∂y ⎥ ⎦ ⎣ ∂x = A11 T x ′ST x + A16 Ty ′ST x + A16 T x ′ST y + A66 Ty ′ST y Kuv

Kvu

K K′c ⎤ K=⎡ s . ⎢Kc 0 ⎥ ⎣ ⎦

∂′ ∂ ∂′ ∂ ∂′ ∂ ∂′ ∂ = ∫ ⎡ A12 + + + A26 A16 A66 ⎤ d  ⎢ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ x y y y x x y x⎥ ⎣ ⎦

The pre-buckling problem of VAT panels is solved by using the equation Ku = F for the displacement vector, u . From the obtained displacements, the stress resultant distributions of the VAT panel can be calculated from Eqs. (2) and (3), which are used in buckling analysis of the VAT panel.

= A12 T x ′ST y + A26 Ty ′ST y + A16 T x ′ST x + A66 Ty ′ST x = Kuv

∂′ ∂ ∂′ ∂ ∂′ ∂ ∂′ ∂ + + + Kvv = ∫ ⎡ A22 A26 A26 A66 ⎤ d  ⎢ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ y y x y y x x x⎥ ⎦ ⎣ = A22 Ty ′ST y + A26 T x ′ST y + A26 Ty ′ST x + A66 T x ′ST x

3.0.2. Buckling analysis

∂′ 1 ∂′ Kuw = ∫ ⎡ ⎛ A12 + A26 ⎞ ⎤ d  ⎢ R ∂x ⎥ ∂y ⎠⎦ ⎣ ⎝ ⎜



By applying the Rayleigh–Ritz method to Eq. (4), that governs the buckling behavior of the VAT panel, the system of equilibrium equations can be written as an eigenvalue problem:

1 = ⎜⎛A12 T x ′S + A26 Ty ′S⎟⎞ R ⎝ ⎠

[K − Pcr K G] u = 0,

Kwu = Kuw ⎜



∂′ ∂ ∂′ ¯ ∂ ∂′ ∂ K G = ∫ ⎡ N¯x Ny + + 2 N¯xy ⎤ d  ⎢ ∂y ∂y ∂x ∂y ⎥ ⎣ ∂x ∂x ⎦ = N¯x T x ′ST x + N¯ y Ty ′ST y + 2N¯xy T x ′ST y,

1 = ⎜⎛A22 Ty ′S + A26 T x ′S⎟⎞ R ⎝ ⎠ = Kvw = ∫ A22 / R2d 

and where, K , K G represents the structural and geometric stiffness matrices, respectively and N¯x , N¯ y and N¯xy represents the pre-buckling stress resultants obtained when a uniform compression load is applied along the edge of the panel.

∂2 ′ ∂2 ∂2 ′ ∂2 ∂2 ′ ∂2 + ∫ ⎡ 2 D11 2 + 2 D12 2 + 2 D16 2 ⎢ ∂x ∂x ∂y ∂x ∂x ∂y ∂x ⎣ ∂2 ′ ∂2 ∂2 ′ ∂2 ∂2 ′ ∂2 + 2 D12 2 + 2 D22 2 + 2 D26 2 ∂x ∂y ∂y ∂y ∂x ∂y ∂y 2′

+2

(22)

where,

∂′ 1 ∂′ Kvw = ∫ ⎡ ⎛ A22 + A26 ⎞ ⎤ d  ⎥ ⎢ R ∂y ∂x ⎠⎦ ⎣ ⎝

Kwv Kww

(21)

∂2

2′

∂2

2′

3.0.3. Boundary conditions

∂2

∂ ∂ ∂ ⎤ d + 2 2 D26 +4 D16 D66 ∂x 2 ∂x ∂y ∂y ∂x ∂y ∂x ∂y ∂x ∂y ⎥ ⎦

In-order to solve the equilibrium equations, corresponding plate boundary conditions have to be specified along the edges. The designation of boundary regions are shown in Fig. 4. The shell ends (i.e. at −a a x ∈ { 2 , 2 } ) are designated as Φx1 and Φx2 , and its sides (i.e. at

= A22 / R2 + D11 T xx ′ST xx + D12 Tyy ′ST xx + 2D16 T xy ′ST xx + D12 T xx ′ST yy + D22 Tyy ′ST yy + 2D26 T xy ′ST yy + 2D16 T xx ′ST xy + 2D26 Tyy ′ST xy + 4D66 T xy ′ST xy,

−b

b

y ∈ { 2 , 2 } ) are designated as Φy1 and Φy2 . Three different types of plate boundary conditions considered in this work, they are: (i) all edges simply supported (SSSS): w (Φxi ) = w (Φyi ) = 0 , (ii) simply supported

(20)

where, ( ) ′ represents the transpose of a given matrix. The pre-buckling analysis of a curved VAT panel is carried out by applying the in-plane boundary conditions: u x = vx = 0atx = −a/2 , vx = 0atx = a/2 and a uniform compression is applied along the edge x = a/2 as shown in Fig. 3. The uniform compression is applied by connecting all nodes on the edge to a control node and an equivalent compression force (fc ) is applied on the control node point (see Fig. 3). The uniform displacement constraint is implemented in GDIQM by

sides (Φy1 and Φy2) and clamped ends (Φx1 and Φx2) (SCSC):

w (Φxi ) = w (Φyi ) = 0 (iii) clamped on all the edges (CCCC): dw (Φyi ) dy

w (Φxi )

w (Φyi )

dw (Φxi ) dx dw (Φxi ) dx

= 0, = 0,

= = 0 , i = 1, 2 . The above boundary condi= 0, tions are implemented by writing them in the GDIQM form. Eq. (22) is solved along with the plate boundary conditions to obtain the critical buckling load Pcr and their corresponding mode shapes. 305

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

following quadratic eigenvalue problems:

[K − (α 0 + 0.5α1) Pcr K G]{a1} − 0.25Mθ 2 {a1} = 0,

(28)

[K − (α 0 − 0.5α1) Pcr K G]{b1} − 0.25Mθ 2 {b1} = 0.

(29)

The solution of Eqs. (28)–(29) determines two boundaries of the first dynamic instability region (principle parametric resonance). 5. Results and discussion In this section, the dynamic instability analysis of a symmetric curved VAT panel subjected to a periodic in-plane compression load is presented. An in-house MATLAB program was developed to analyze the curved VAT panel using GDIQM. For numerical study, a symmetric curved VAT panel [ϕ ± 〈T0 |T1 〉]2s of dimensions a = b = 0.3 m and thickness h = 1.048 × 10−3 m is considered. The material properties of the lamina are given as: E1 = 163 GPa, E2 = 10 GPa, G12 = 5 GPa, ν12 = 0.3, ρ = 1480 kg/m3 . Simply supported boundary conditions has been used for all analyses unless otherwise specified. Firstly, the accuracy of the GDIQM is verified by comparison with results published in the literature and with ABAQUS simulations. For ABAQUS, a subroutine was developed to generate four noded shell elements (S4R) of curved VAT panel, each having independent fiber orientations. To study the convergence of GDIQM, the normalized fundamental ωa2 ρ / E2 ) of the curved VAT panels ([90 ± 〈0|0〉]2s and frequency (Ω˜ =

Fig. 4. Representation of boundaries of a curved panel with VAT laminate.

4. Dynamic stability of curved VAT panel Dynamic stability analysis is performed on a curved VAT panel subjected to a periodic axial compression load of form P (t ) = α 0 Pcr + α1 Pcr cos(θt ) . Here, Pcr is the critical buckling load of the panel, α 0 , α1 are the static and dynamic load parameters, respectively and θ is the forcing frequency. The equations of motion of a symmetric curved VAT panel under periodic axial compression load is expressed as:

Mu ¨ + Ku − (α 0 Pcr + α1 Pcr cos(θt )) K G u = 0,

h

[90 ± 〈0|30〉]2s ) with b/ R = 0.2 for different number of grid points (n x × n y ) is shown in Fig. 5. FEM results were computed using ABAQUS for different mesh sizes (see Fig. 5) to obtain convergence. From Fig. 5 it is observed that a grid size of n x = n y = 21 is sufficient to obtain accurate results and therefore, the same grid size has been used for all analyses in this work. Also, it is observed that GDIQM requires fewer grid points when compared to FEM (ABAQUS) to obtain converged results. This study shows the fast convergence and less computational effort required by GDIQM than FEM. Modal analysis of a curved VAT panel was undertaken using GDIQM and the results are now compared with ABAQUS simulations. The analysis is performed on a curved VAT panel [90 ± 〈0|T1 〉]2s for different radius of curvature and the corresponding normalized fundamental natural frequencies along with the ABAQUS simulation results are presented in Table 1. The first four vibration mode shapes of a curved VAT panel [90 ± 〈0|30〉]2s with b/ R = 0.2 obtained from the present analysis are compared with the ABAQUS results in Table 2. As shown in Table 1 and 2, normalized natural frequencies and mode shapes obtained using the GDIQM matches closely with the ABAQUS results.

(23)

where, M represents the mass matrix of the VAT panel and is obtained from the kinetic energy expression of the panel, which is given as:

KE =

1 2

∫ ρ [h (u˙ 2 + v˙ 2 + w˙ 2)] d = 12 u˙ T M u˙ ,

(24)

here, ρ is the mass density and h is the thickness of the panel. From Eq. (24), the mass matrix for the panel can be written in GDIQM form as:

⎡ hS 0 0 ⎤ M = ρ ⎢ 0 hS 0 ⎥ ⎢ ⎦ ⎣ 0 0 hS ⎥

(25)

In Eq. (23), for certain combinations of α 0 , α1 and θ , the response of the panel can become unbounded due to parametric resonance and is referred to as dynamic instability. The region in the parametric space, where the dynamic instability occurs is known as the dynamic instability region (DIR). In this work, Bolotin's approach [12] has been used to determine the stability boundary of Eq. (23). According to Bolotin's theory, the solutions of Eq. (23) are periodic with period T1 = 2π / θ or T2 = 4π / θ on the stability boundary. Since, the principle dynamic instability region is the most dominant among all other instability regions and has the boundary with the solutions of period 4π / θ , the solution of Eq. (23) can be expanded as: ∞

u=

∑ k = 1,3,5

⎛ ⎞⎞ ⎛ ⎛ ⎞ ⎜a k sin ⎜kθt /2⎟ + bk cos ⎜kθt /2⎟ ⎟, ⎝ ⎠ ⎝ ⎠⎠ ⎝

(26)

where, the coefficients a k , bk are unknowns and infinite in number. Therefore, a first-order approximation is often used in the literature [14,37,19] to predict the stability boundary. In the first-order Bolotin's approximation, only the first two terms in Eq. (26) are considered, therefore we get:

u = a1 sin(θt /2) + b1 cos(θt /2).

(27) Fig. 5. Numerical convergence study of normalized fundamental frequency (Ω˜ ) computed using GDIQM and FEM (ABAQUS) for the curved VAT panel with b/ R = 0.2 .

Upon substituting Eq. (27) into Eq. (23) we obtain the residual equations. It should be noted that the coefficients of sin(θt /2) and cos(θt /2) must vanish in the residual equations. This leads to the 306

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

Table 1 Comparison of normalized fundamental natural frequency (Ω˜ ) of a curved VAT panel [90 ± 〈0|T1 〉]2s for different b/ R ratios. T1 b/R

0 0.05 0.075 0.1 0.2



60°

30°

Present

Abaqus

Present

Abaqus

Present

Abaqus

2.0279 2.3779 2.7493 3.1951 5.3125

2.0310 2.3802 2.7536 3.2013 5.3287

2.3176 2.9177 3.5057 4.1833 7.2526

2.3146 2.9184 3.5113 4.1928 7.2799

2.5294 3.5526 4.5065 5.5681 7.3138

2.5229 3.5526 4.5130 5.5791 7.3228

Next, the buckling load of a VAT panel for a range of curvatures (b/ R = 0.0, 0.1, 0.2) and fiber angle orientations are determined using GDIQM and is presented in Fig. 6 along with the results published in the literature [38] and ABAQUS simulation results. Fig. 6 clearly illustrates a very good agreement between the present model and the results published in the literature and ABAQUS results. The comparison of first four buckling mode shapes along with the normalized critical buckling load values (K cr = Pcr a2 /(E1 bh3)) of a curved VAT panel [90 ± 〈0|30〉]2s with b/ R = 0.2 from the present model with ABAQUS is shown in Table 3. We also study, the variation of normalized natural frequency with respect to the static in-plane load parameter (α 0 ) for a curved VAT panel [90 ± 〈0|30〉2s ] with b/ R = 0.1. Here, the VAT panel is subjected to a uniform axial compression load with the static load parameter, α 0 varying from 0.0 to 1.0 and the corresponding normalized natural frequencies are plotted in Fig. 7. The results are then compared with ABAQUS and found to be in good agreement with each other.

Fig. 6. Variation of buckling load of a curved VAT panel [0 ± 〈0|T1 〉]2s with respect to fiber angle orientation T1 .

T1 = 60° and then decreases as shown in Fig. 8. This change in the parametric frequency and critical buckling load occurs as a consequence of the change in membrane and bending stiffness of the VAT panel due to the change in fiber angle orientations. Now, the influence of radius of curvature (b/ R) on DIR of a curved VAT panel is studied. The analysis is carried out on a VAT panel [90 ± 〈0|30〉]2s for four different b/ R ratios (0.0, 0.05, 0.075, 0.1, 0.2) and the results are presented in Fig. 9. The ratio b/ R = 0 in Fig. 9 corresponds to a flat laminate. In Fig. 9, for the VAT panel when b/ R = 0 , the onset of DIR starts at Ω˜ = 4.97 and extends up to K cr = 3.92. In the case of b/ R = 0.2, the onset of DIR shifts to Ω˜ = 15.29 and extends up to K cr = 13.03, thus increasing the dynamic stability region of the VAT panel. Also, the width of the DIR shrinks with increase in the b/ R ratio of the panel. This study clearly illustrates the increase in dynamic stability of the curved VAT panel with increase in the b/ R ratio.

5.1. Dynamic instability of curved VAT panel In this section, the effect of fiber angle variation on dynamic instability region (DIR) of a curved VAT panel is studied. The analysis is performed on a VAT panel [90 ± 〈0|T1 〉]2s for different values of T1 ranging from 0° to 90° in increments of 15°. The radius of curvature of the panel is taken as b/ R = 0.2 and the static load parameter is chosen as α 0 = 0.0 for the analysis. The results are shown in Fig. 8 in the plane of normalized frequency (Ω˜ ) versus normalized buckling load (K cr = Pcr a2 /(E1 bh3)) . From Fig. 8 it is clearly seen that a change in the fiber angle configuration influences the DIR of a VAT panel. The onset of DIR shifts towards higher frequency when T1 increases, reaches a maximum when T1 = 45° and then decreases. Likewise, the extension of DIR towards the K cr value increases with T1, reaches a maximum when

5.2. Variation of buckling load and dynamic instability index of a curved VAT panel In this section, the variation of critical buckling load and dynamic instability index of a curved VAT panel with respect to fiber angle orientation is studied. The analysis is carried out on a VAT panel [90 ± 〈T0 |T1 〉]2s with radius of curvature b/ R = 0.2. The variation of K cr with respect to the fiber angle T1 for different values of T0 is shown in Fig. 10(a). Each curve is generated by varying T1 from 0° to 90° with an increment of 5° by considering T0 as a constant value. Variation of K cr

Table 2 Comparison of free vibration mode shapes of a curved VAT panel [90 ± 〈0|30〉]2s .

307

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

Table 3 Comparison of buckling mode shapes of a curved VAT panel [90 ± 〈0|30〉]2s .

Fig. 7. Variation of normalized fundamental natural frequency (Ω˜ ) with respect to static in-plane load parameter, α 0 for a curved VAT panel [90 ± 〈0|30〉]2s .

Fig. 9. Dynamic instability regions of a curved VAT panel [90 ± 〈0|45〉]2s for different b/ R ratios.

majority of the applied compression load is redistributed from the center towards the supported edges by the fiber orientation paths. Likewise, in the case of curved VAT panels, the benign load redistribution enables VAT layups to attain high buckling resistance. The load redistribution in curved VAT panels is governed by the fiber angle orientations and the radius of curvature. From Fig. 10(a) it is observed that the redistribution of applied compression load from the buckling regions to the supported edges is attained maximum for the fiber angle orientation T0 = 15° and then decreases with increase in T0 from 15° to 90°. Moreover, for a constant T0 , the normalized critical buckling load first increases with increase in the fiber angle T1 and then decreases. In Fig. 10(a), among all the straight fiber laminate configurations, [± 60]2s shows maximum normalized critical buckling load value of 7.6. In the case of VAT panels, the highest normalized K cr value is 11.48 for the fiber angle orientation T0 = 15° and T1 = 60°. Thus, the VAT panel ([90 ± 〈15|60〉]2s ) exhibits 51% improvement in buckling performance when compared to the optimal straight fiber laminate ([± 60]2s ). Also, VAT panels with the fiber angle orientations between T0 = 0°to45° shows higher critical buckling load values for a range of T1 = 45°to90°, when compared to straight fiber laminates. Next, the dynamic instability of a curved VAT panel for different fiber angle orientations is studied. Since, dynamic instability of a VAT panel is closely related to its fundamental frequency and critical buckling load, a non-dimensional parameter, dynamic instability index (DII ) is used to characterize the dynamic stability [21,19] and is given as:

Fig. 8. Dynamic instability regions of a curved VAT panel [90 ± 〈0|T1 〉]2s for different fiber orientations, T1 .

for a straight fiber laminate ([±θ]2s ) is also shown in Fig. 10(a). From the results, it is clearly observed that the K cr value of a VAT panel changes significantly with the fiber angle orientations (T0 andT1) . This phenomenon arises due to the variation in the membrane (A) and bending (D) stiffness of the VAT panel with the linearly varying fiber paths. As shown in previous works [3,4], on flat VAT laminate, the 308

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

Fig. 10. Variation of normalized critical buckling load (K cr ) and DII for the SSSS curved VAT panel [90 ± 〈T0 |T1 〉]2s with b/ R = 0.2 for different fiber orientations, T0 and T1.

DII =

DIO , ωK cr

CCCC (clamped on all the edges) boundary conditions and the results are shown in Fig. 11 and 12, respectively. In the SCSC case, the optimal VAT panel with fiber angles T0 = 15° and T1 = 65° shows 44% higher normalized critical buckling load than the optimal straight fiber laminate with layup [± 60]2s . Also, the minimum value of DII observed in the case of a VAT panel is 0.01835 for the fiber angles T0 = 0° and T1 = 70°, which is 20% less than the optimal straight fiber laminate [± 30]2s . Whereas, in the case of CCCC boundary condition, a curved VAT panel with configuration [90 ± 〈15|60〉]2s shows 43.79% higher normalized critical buckling load value and a VAT panel of [90 ± 〈0|65〉]2s shows 9.22% decrease in DII value when compared to the optimal straight fiber laminate [± 45]2s . Furthermore, Fig. 10, 11 and 12 collectively show that many VAT panel configurations with different fiber angle orientations exhibits higher buckling resistance and dynamic stability than straight fiber laminates, thus providing additional freedom in tailoring the stiffness of the laminate. The results clearly demonstrate the advantages and the potential of applying variable angle tow fiber paths to improve the buckling and dynamic performance of a curved VAT panel. As an additional consideration, DIR of a VAT panel with fiber angle configuration [90 ± 〈0|30〉]2s and b/ R = 0.2 for different boundary conditions (SSSS, SCSC and CCCC) is shown in Fig. 13. The static load parameter α 0 is chosen to be zero for the analysis. From Fig. 13 it is observed that from SSSS to CCCC the onset of DIR shifts towards higher Ω˜ and extends towards larger K cr values. This behavior occurs because of the increase in both natural frequency and critical buckling load of the VAT panel

(30)

where, DIO is the width of the dynamic instability region evaluated at a specific value of α1. A higher value of DII indicates that the VAT panel is more dynamically unstable. The variation of DII with respect to fiber angle T1 is shown in Fig. 10(b). The results are generated by varying T1 from 0° to 90° with an increment of 5° for a given value of T0 . The static and dynamic load parameters are chosen as α 0 = 0 and α1 = 0.3 for the analysis. The results clearly show that the variation in fiber angle configuration of the VAT panel affects the DII value considerably. From Fig. 10(b) it is observed that many of the VAT panel configurations from T1 = 40° to T1 = 90° show lower DII values when compare to straight fiber laminates. The minimum value of DII observed for a straight fiber laminate is 0.02155 for the layup [± 30]2s . However, in the case of a VAT panel, the minimum value of DII is found to be 0.0188 for the fiber angle orientations T0 = 0 and T1 = 65°, and this value is 14% less than the minimum value obtained for optimal straight fiber laminate configuration. This substantial decrease in DII value occurs because of the increase in critical buckling load and the natural frequency by the linearly varying fiber orientation paths of the VAT panel. This increase in critical buckling and natural frequency reduces the DII value, thus makes the VAT panel dynamically more stable. A similar analysis is carried out on a curved VAT panel [90 ± 〈T0 |T1 〉]2s with SCSC (clamped sides and simply supported ends),

Fig. 11. Variation of normalized critical buckling load and DII for the SCSC curved VAT panel [90 ± 〈T0 |T1 〉]2s with b/ R = 0.2 for different fiber orientations, T0 and T1. 309

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

Fig. 12. Variation of normalized critical buckling load and DII for the CCCC curved VAT panel [90 ± 〈T0 |T1 〉]2s with b/ R = 0.2 for different fiber orientations, T0 and T1.

Fig. 15. Variation of K cr of a curved VAT panel [90 ± 〈0|T1 〉]2s with respect to T1 for different values of a/ b ratios.

Fig. 13. Dynamic instability region of a curved VAT panel [90 ± 〈0|30〉]2s with b/ R = 0.2 for different boundary conditions (SSSS, SCSC and CCCC).

particular b/ R ratio. From the results, it can be deduced that with an increase in b/ R ratio (i.e decrease in radius of curvature), the critical buckling load increases for all fiber angle orientations of a VAT panel. For a particular b/ R ratio, K cr value increases when the fiber angle T1 increases from 0° to 70° and then decreases. This behavior happens

under clamped boundary conditions. The variation of K cr and DII of a curved VAT panel with respect to the radius of curvature is shown in Fig. 14(a). The analysis is carried out on a VAT panel [90 ± 〈0|T1 〉]2s for different values of T1. In Fig. 14(a) each curve is generated by varying the fiber angle T1 from 0° to 90° for a

Fig. 14. Variation of K cr and DII for a SSSS VAT panel [90 ± 〈0|T1 〉]2s with respect to T1 different values of b/ R ratios. 310

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

Fig. 16. Variation of DIR and DII of a curved VAT panel [90 ± 〈0|T1 〉]2s with respect to T1 for different values of a/ b ratios.

because of the redistribution of the in-plane stresses by the fiber paths. This redistribution of stresses from the center of the panel towards the supported edges increases until the fiber angle T1 reaches 70° and then decreases. The variation of DII with respect to fiber angle T1 is also shown in Fig. 14(b). In this case, the value of DII decreases with an increase in b/ R ratio (decrease in radius of curvature) of the panel and clearly indicates that the buckling phenomenon governs the dynamic stability of a VAT panel. Finally, the variation of K cr with respect to the fiber angle T1, for different aspect ratios (a/ b) is studied. The analysis is done on a VAT panel [90 ± 〈0|T1 〉]2s for different values of a/ b ratio (0.5, 1.0, 2.0, 3.0) by varying fiber angle T1 from 0° to 90°. The analysis is carried out by fixing the width (b = 100 × 10−3m) and varying the length of the panel. The radius of curvature of the panel is taken as (b/ R = 0.2) and the results are shown in Fig. 15. From the results, it is observed that with the increase in a/ b ratio, normalized critical buckling load of the VAT panel increases significantly for all fiber angle configurations. For a particular value of a/ b ratio the value of K cr increases when the fiber angle T1 increases from 0° to 70° and then decreases. DIRs of a VAT panel [90 ± 〈0|45〉]2s for different values of a/ b ratio is shown in Fig. 16(a). It can be seen from the figure that with an increase in a/ b ratio DIR shifts to higher parametric frequency and extends to larger K cr value, indicating an increase in critical buckling load and natural frequency of the VAT panel. Also, the width of DIR shrinks with the increase in a/ b ratio, illustrating the increase in dynamic stability of the panel. The variation of DII with respect to the fiber angle T1 for different a/ b ratios is shown in Fig. 16(b). As expected, the DII value of a VAT panel decreases with increase in a/ b ratio, because of the increase in critical buckling load and natural frequency of the panel. From all results presented it is clearly observed that the curved panel with VAT laminate for some fiber angle configurations shows better buckling and dynamic performance when compared to conventional optimal straight fiber laminate. Moreover, the steering of fiber orientation over the curved lamina provides an added advantage to the designer by allowing more scope for tailoring K cr , DII and stiffness of the VAT laminate simultaneously.

Initially, the effect of tow steering on buckling and dynamic instability behavior of a curved VAT panel for three different types of boundary conditions was studied. It is observed that the fiber angle variation has a significant effect on critical buckling load as well as dynamic instability of the curved VAT panel. For some fiber angle configurations, curved VAT panel show much higher critical buckling load and dynamic stability than straight fiber laminates. This improvement in performance of a curved VAT panel is because of the redistribution of the applied load from the center of the panel towards the supported edges by the curvilinear fiber paths. In addition, the influence of radius of curvature on the behavior of a curved VAT panel has been presented. It is observed that the critical buckling load and dynamic stability increases with the decrease in the radius of curvature of a curved VAT panel irrespective of the fiber angle configurations. Also, with the increase in the aspect ratio, buckling and dynamic performance of a curved VAT panel improves. From this study, the benefits of fiber steering in improving the buckling and dynamic behavior of a curved VAT panel has been demonstrated. Acknowledgment C.P. Vyasarayani gratefully acknowledges the Department of Science and Technology India for funding this research through the Inspire fellowship (DST/INSPIRE/04/2014/000972). Paul M. Weaver would like to acknowledge the support of Science Foundation Ireland for the award of a Research Professor grant (Varicomp: 15/RP/2773). References [1] M.W. Hyer, H. Lee, The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes, Compos. Struct. 18 (3) (1991) 239–261. [2] Z. Gurdal, R. Olmedo, In-plane response of laminates with spatially varying fiber orientations-variable stiffness concept, AIAA J. 31 (4) (1993) 751–758. [3] Z. Gurdal, B.F. Tatting, C. Wu, Variable stiffness composite panels: effects of stiffness variation on the in-plane and buckling response, Compos. Part A: Appl. Sci. Manuf. 39 (5) (2008) 911–922. [4] G. Raju, Z. Wu, B.C. Kim, P.M. Weaver, Prebuckling and buckling analysis of variable angle tow plates with general boundary conditions, Compos. Struct. 94 (9) (2012) 2961–2970. [5] Z. Wu, P.M. Weaver, G. Raju, B.C. Kim, Buckling analysis and optimisation of variable angle tow composite plates, Thin-Walled Struct. 60 (2012) 163–172. [6] Z. Wu, G. Raju, P.M. Weaver, Postbuckling analysis of variable angle tow composite plates, Int. J. Solids Struct. 50 (10) (2013) 1770–1780. [7] G. Raju, Z. Wu, P.M. Weaver, Buckling and postbuckling of variable angle tow composite plates under in-plane shear loading, Int. J. Solids Struct. 58 (2015) 270–287. [8] B.H. Coburn, P.M. Weaver, Buckling analysis, design and optimisation of variable-

6. Conclusion In this work, the dynamic stability behavior of a curved panel with VAT laminate subjected to periodic axial compression load has been studied. Donnell's shallow shell theory was used to model the curved VAT panel and solved using the generalized differential integral quadrature method. Bolotin's first-order approximation was used to determine the dynamic instability region of a curved VAT panel. 311

Thin-Walled Structures 138 (2019) 302–312

S. Samukham, et al.

(1994) 335–342. [24] S. Sahu, P. Datta, Parametric resonance characteristics of laminated composite doubly curved shells subjected to non-uniform loading, J. Reinf. Plast. Compos. 20 (18) (2001) 1556–1576. [25] T. Ng, K. Lam, J. Reddy, Parametric resonance of a rotating cylindrical shell subjected to periodic axial loads, J. Sound Vib. 214 (3) (1998) 513–529. [26] K. Lam, T. Ng, Dynamic stability of cylindrical shells subjected to conservative periodic axial loads using different shell theories, J. Sound Vib. 207 (4) (1997) 497–520. [27] L.H. Donnell, Stability of thin-walled tubes under torsion, No.NACA-R-479 California Inst. of Tech., Pasadena, USA. [28] G. Raju, Z. Wu, P.M. Weaver, Postbuckling analysis of variable angle tow plates using differential quadrature method, Compos. Struct. 106 (2013) 74–84. [29] S.C. White, G. Raju, P.M. Weaver, Initial post-buckling of variable-stiffness curved panels, J. Mech. Phys. Solids 71 (2014) 132–155. [30] G. Raju, Z. Wu, S. White, P.M. Weaver, Optimal postbuckling design of variable angle tow composite plates, AIAA J. 56 (5) (2018) 2045–2061. [31] R. Bellman, B. Kashef, J. Casti, Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations, J. Comput. Phys. 10 (1) (1972) 40–52. [32] C. Shu, Differential Quadrature and Its Application in Engineering, Springer Science & Business Media, 2012. [33] A. Sherbourne, M. Pandey, Differential quadrature method in the buckling analysis of beams and composite plates, Comput. Struct. 40 (4) (1991) 903–913. [34] X. Wang, M. Tan, Y. Zhou, Buckling analyses of anisotropic plates and isotropic skew plates by the new version differential quadrature method, Thin-Walled Struct. 41 (1) (2003) 15–29. [35] R. Penrose, A generalized inverse for matrices, Mathematical Proceedings of the Cambridge Philosophical Society, 51 Cambridge University Press, 1955, pp. 406–413. [36] H. Zhang, F. Ding, On the kronecker products and their applications, J. Appl. Math. (2013). [37] A.G. Radu, A. Chattopadhyay, Dynamic stability analysis of composite plates including delaminations using a higher order theory and transformation matrix approach, Int. J. Solids Struct. 39 (7) (2002) 1949–1965. [38] S.C. White, Post-buckling of variable-stiffness shell structures, Ph.D. Thesis, University of Bristol.

stiffness sandwich panels, Int. J. Solids Struct. 96 (2016) 217–228. [9] P. Hao, C. Liu, X. Yuan, B. Wang, G. Li, T. Zhu, F. Niu, Buckling optimization of variable-stiffness composite panels based on flow field function, Compos. Struct. 181 (2017) 240–255. [10] B.C. Kim, K. Potter, P.M. Weaver, Continuous tow shearing for manufacturing variable angle tow composites, Compos. Part A: Appl. Sci. Manuf. 43 (8) (2012) 1347–1356. [11] H.-J.L. Dirk, C. Ward, K.D. Potter, The engineering aspects of automated prepreg layup: history, present and future, Compos. Part B: Eng. 43 (3) (2012) 997–1009. [12] V. Bolotin, The Dynamic Stability of Elastic Systems, Holden-Day, San Francisco, 1964. [13] R. Srinivasan, P. Chellapandi, Dynamic stability of rectangular laminated composite plates, Comput. Struct. 24 (2) (1986) 233–238. [14] J. Moorthy, J. Reddy, R. Plaut, Parametric instability of laminated composite plates with transverse shear deformation, Int. J. Solids Struct. 26 (7) (1990) 801–811. [15] A. Chattopadhyay, A.G. Radu, Dynamic instability of composite laminates using a higher order theory, Comput. Struct. 77 (5) (2000) 453–460. [16] P. Dey, M. Singha, Dynamic stability analysis of composite skew plates subjected to periodic in-plane load, Thin-walled Struct. 44 (9) (2006) 937–942. [17] S. Sahu, P. Datta, Dynamic instability of laminated composite rectangular plates subjected to non-uniform harmonic in-plane edge loading, Proc. Inst. Mech. Eng. Part G: J. Aerosp. Eng. 214 (5) (2000) 295–312. [18] G. Cederbaum, Dynamic instability of shear-deformable laminated plates, AIAA J. 29 (11) (1991) 2000–2005. [19] S. Wang, D. Dawe, Dynamic instability of composite laminated rectangular plates and prismatic plate structures, Comput. Methods Appl. Mech. Eng. 191 (17–18) (2002) 1791–1826. [20] X. Chen, G. Nie, Z. Wu, Dynamic instability of variable angle tow composite plates with delamination, Compos. Struct. 187 (2018) 294–307. [21] S. Samukham, G. Raju, C. Vyasarayani, Parametric instabilities of variable angle tow composite laminate under axial compression, Compos. Struct. 166 (2017) 229–238. [22] V. Khalafi, J. Fazilati, Parametric instability behavior of tow steered laminated quadrilateral plates using isogeometric analysis, Thin-Walled Struct. 133 (2018) 96–105. [23] M. Ganapathi, T. Varadan, V. Balamurugan, Dynamic instability of laminated composite curved panels using finite element method, Comput. Struct. 53 (2)

312