- Email: [email protected]
Dynamic stability analysis of composite laminated cylindrical panels via the mesh-free kp-Ritz method
ARTICLE IN PRESS
International Journal of Mechanical Sciences 49 (2007) 1156–1165 www.elsevier.com/locate/ijmecsci
Dynamic stability analysis of composite laminated cylindrical panels via the mesh-free kp-Ritz method K.M. Liewa,, Y.Y. Leea, T.Y. Ngb, X. Zhaob a
b
Department of Building and Construction, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong School of Mechanical and Production Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore Received 14 June 2006; received in revised form 19 February 2007; accepted 19 February 2007 Available online 24 February 2007
Abstract In this paper, the dynamic stability analysis of thin, laminated cylindrical panels under static and periodic axial forces is presented by using the mesh-free kp-Ritz method. The mesh-free kernel particle estimate is employed to approximate the 2D transverse displacement field. A system of Mathieu–Hill equations is obtained through the application of the Ritz minimization procedure to the energy expressions. The principal instability regions are then analyzed via Bolotin’s first approximation. Effects of lamination schemes such as number of layers, ply-angle and boundary conditions on the instability regions are examined in detail. r 2007 Elsevier Ltd. All rights reserved. Keywords: Dynamic stability; Parametric resonance; Cylindrical panel; Composite laminate; Mesh-free kp-Ritz method; Bolotin’s first approximation
1. Introduction The dynamic stability or phenomenon of parametric resonance of a structure under periodic loads has attracted an enormous research interest. The earlier work include the dynamic stability analysis of elastic system by Bolotin [1], the parametric instability of circular cylindrical shells by Vijayaraghavan and Evan-Iwanowski [2], and nonlinear elastic buckling and parametric excitation of a cylinder under axial loads by Yao [3]. Birman [4] and Srinivasan and Chellapandi [5] studied the dynamic stability of laminated rectangular plates. The parametric instability of plates with transverse shear deformation was investigated by Moorthy and Reddy [6] and Bert and Birman [7]. Bert and Birman [8] extended Yao’s [3] approach to the parametric instability of thick orthotropic shells using higher-order theory. Koval [9] analyzed the effects of longitudinal resonance on the parametric stability of an axially excited cylindrical shell. Argento and Scott [10,11] employed a perturbation technique to examine the dynamic stability of layered anisotropic circular cylindrical shells Corresponding author. Tel.: +852 3442 65816; fax: +852 2788 7612.
E-mail address: [email protected] (K.M. Liew). 0020-7403/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2007.02.005
under axial loading. More recently, the dynamic stability of cross-ply laminated composite cylindrical shells under combined static and periodical axial forces was investigated by Ng et al. [12] and Lam and Ng [13]. The dynamic instability of laminated composite curved panels was presented by Ganapathi et al. [14] using finite element method. Ng et al. [15] investigated the dynamic stability of a simply supported cylindrical panel using generalizes Donnell’s shell theory. The Ritz or Rayleigh–Ritz method is a proven approximate technique in computational mechanics, and notable works include those of Kitipornchai et al. [16], Liew et al. [17,18], Cheung and Zhou [19,20], Zeng and Bert [21], Su and Xiang [22], Lim and Liew [23], and Liew and Feng [24]. In this paper, for the dynamic stability analysis of composite laminated cylindrical shell panels under combined static and periodic axial loads, an energy formulation is first described, and a system of Mathieu–Hill equations is obtained via the Ritz minimization procedure. The parametric resonance responses are then analyzed using Bolotin’s [1] method. In conventional Ritz method, it is difficult to choose appropriate functions to satisfy some complicate boundary conditions, and the eigen equations need to be reevaluated for different boundary condition
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
cases. In the present kp-Ritz method, a standard weight function is employed to express the interior field, and the boundary conditions are enforced by penalty method. Therefore, the disadvantages in conventional Ritz method are avoided. For the assumed functions in the present work, the mesh-free kernel particle approximate is utilized to describe the 2D transverse displacement field. The commonly used trigonometric or hyperbolic functions in conjunction with the Ritz procedure are not suitable for general laminates as erroneously vanish. The present formulation ensures that the bending–extension coupling stiffnesses in general laminates are retained. The present method is verified by comparing with those results in open literature. The effects of boundary conditions and lamination schemes on the instability regions are also investigated. 2. Theoretical formulation
1157
load is given by N a ¼ N 0 þ N s cos Pt,
(1)
where P is the frequency of excitation in radians per unit time. The kinetic energy for the cylindrical panels can be expressed as Z L Z y0 1 ½u_ 2 þ v_2 þ w_ 2 R dy dx, (2) T ¼ rh 2 0 0 where the three terms are due to contributions from the linear velocities in the x, y, and z directions, respectively. The strain energy U a due to the axial loading can be written as " 2 2 # Z Z 1 L y0 qu 2 qv qw Ua ¼ Na þ þ R dy dx. 2 0 0 qx qx qx (3)
2.1. Energy formulation for cylindrical panels The cylindrical panel considered is shown in Fig. 1a, where a coordinate system ðx; y; zÞ is fixed on the middle surface of the panel. This panel is considered to be thin and of length L, radius R and thickness h, bounded along its edges by the lines x ¼ 0, x ¼ L, y ¼ 0 and y ¼ y0 . The displacements of the panel in the x, y, and z directions are denoted by u, v and w, respectively. The pulsating axial
s
L
w
h
v θο
b
u
R θ
T ¼ fe1 e2 g k1 k2 2tg,
(5)
where the middle surface strains, e1 , e2 and g, and the middle surface curvatures, k1 , k2 and t, are defined according to Love’s thin shell theory qu 1 qv qv 1 qu ; e2 ¼ þw ; g¼ þ e1 ¼ qx R qy qx R qy
and [S] is given by 2 A11 A12 A16 6 6 A12 A22 A26 6 6 A16 A26 A66 6 ½S ¼ 6 6 B11 B12 B16 6 6 B12 B22 B26 4 B16 B26 B66
v
h2
h4 h5
h
B11
B12
B16
B12 B16
B22 B26
B26 B66
D11
D12
D12 D16
D22 D26
ð6Þ
3
7 7 7 7 7 7, D16 7 7 D26 7 5 D66
(7)
where the extensional ðAij Þ, coupling ðBij Þ and bending ðDij Þ stiffnesses, are defined as Z
Fig. 1. (a) Coordinate system of the cylindrical panel; (b) cross-sectional view of the laminated cylindrical panel.
(4)
where T and [S] are the strain vector and stiffness matrix, respectively, and T can be defined as
q2 w 1 q2 w qv k1 ¼ 2 ; k2 ¼ 2 , R q2 y qy q x 1 q2 w qv t¼ R qxqy qx
w
h1
The strain energy of the shell is given by Z Z 1 L y0 T U ¼ ½SR dy dx, 2 0 0
h=2
ðAij ; Bij ; Dij Þ ¼ h=2
Qij ð1; z; z2 Þ dz.
(8)
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1158
For a panel composed of different layers of materials, the stiffnesses are expressed as Nl X ¯ kij ðhk hkþ1 Þ, Q Aij ¼ k¼1
Bij ¼
l 1X ¯ k ðh2 h2kþ1 Þ, Q 2 k¼1 ij k
Dij ¼
Nl 1X ¯ k ðh3 h3 Þ, Q kþ1 3 k¼1 ij k
N
cI ðxÞ ¼ Cðx; x xI ÞFa ðx xI Þ,
ð9Þ
Cðx; x xI Þ ¼ H T ðx xI ÞbðxÞ,
cos2 a sin2 a
where a is the angular orientation of the fibres and [Q] is the reduced stiffness matrix defined as 2 3 Q11 Q12 0 6 0 7 (12) ½Q ¼ 4 Q12 Q22 5. 0 0 Q66 The material constants in the reduced stiffness matrix [Q] are given as E 11 n12 E 22 Q11 ¼ ; Q12 ¼ , 1 n12 n21 1 n12 n21 E 22 Q22 ¼ ; Q66 ¼ G 12 , ð13Þ 1 n12 n21 where E 11 and E 22 are the elastic moduli in the principle material coordinates, G 12 is the shear modulus and n12 and n21 are Poisson’s ratios. The total energy functional of the panel is thus Gt ¼ T U a U . (14)
bT ðxÞ ¼ ½b0 ðx; yÞ; b1 ðx; yÞ; b2 ðx; yÞ; b3 ðx; yÞ; b4 ðx; yÞ; b5 ðx; yÞ, H T ðx xI Þ ¼ ½1; x xI ; y yI ; ðx xI Þðy yI Þ, ðx xI Þ2 ; ðy yI Þ2 .
2.2. Two-dimensional kernel particle shape functions The discrete displacement approximations for shell panels can be expressed as
cI ðxÞ ¼ bT ðxÞHðx xI ÞFa ðx xI Þ. cI ðxÞ ¼ bT ðxÞB I ðx xI Þ,
vðx; yÞ ¼
where bðxÞ ¼ M 1 ðxÞHð0Þ,
ð21Þ
BI ðx xI Þ ¼ Hðx xI ÞFa ðx xI Þ
ð22Þ
and the moment matrix M is a function of x while Hð0Þ is a constant vector. The expressions of M and Hð0Þ are given by MðxÞ ¼
NP X
Hðx xI ÞH T ðx xÞFa ðx xI Þ,
NP X I¼1
cI ðx; yÞwI eiot ,
ð23Þ
I¼1
Hð0Þ ¼ ½1; 0; 0; 0; 0; 0T .
ð24Þ
The shape function can therefore be expressed as cI ðxÞ ¼ H T ð0ÞM 1 ðxÞHðx xI ÞFa ðx xI Þ.
(25)
For this thin shell problem, the first and second derivatives of shape functions need to be determined. Eq. (21) can be rewritten as (26)
The vector bðxÞ can be determined using LU decomposition of the matrix MðxÞ followed by backward substitution. The derivatives of bðxÞ can be obtained similarly. Taking the derivative of Eq. (26), we have
MðxÞb;x ðxÞ ¼ H ;x ð0Þ M ;x ðxÞbðxÞ.
cI ðx; yÞvI eiot ,
I¼1
wðx; yÞ ¼
(20)
(27)
which can be rearranged as
I¼1 NP X
(19)
Eq. (16) can be rewritten as
M ;x ðxÞbðxÞ þ MðxÞb;x ðxÞ ¼ H ;x ð0Þ
cI ðx; yÞuI eiot ,
ð18Þ
H is a vector of quadratic basis and bi ðx; yÞ’s are functions of x and y which are to be determined. Thus, the shape function can be written as
MðxÞbðxÞ ¼ Hð0Þ.
NP X
(17)
where
where [T] is the transformation matrix for the principle material coordinates and the panel’s coordinates, and is defined as 2 3 sin2 a 2 cos a sin a cos2 a 6 7 ½T ¼ 4 sin2 a (11) cos2 a 2 cos a sin a 5,
uðx; yÞ ¼
(16)
where Cðx; x xI Þ is the correction function and Fa ðx xI Þ is the kernel function. The correction function Cðx; x xI Þ is described as
where hk and hkþ1 denote the distances from the panel reference surface to the outer and inner surfaces of the kth layer, respectively, as shown in Fig. 1b. N l denotes the total ¯ kij is the number of layers in the laminated panel and Q transformed reduced stiffness matrix for the kth layer defined as ¯ ¼ ½T½Q½TT , ½Q (10)
cos a sin a cos a sin a
where NP is the total number of particles and cI ðx; yÞ’s are the shape functions of displacements u, v, and w. The 2D shape function is expressed as (see [25])
ð15Þ
(28)
It can be seen that the first derivative of bðxÞ can be formulated using the LU decomposition procedure again. The second derivative of bðxÞ can be determined by taking derivative of Eq. (28) and using the same LU
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
decomposition procedure. The first derivative of the shape function can be obtained by taking the derivative of Eq. (20): cI;x ðxÞ ¼
bT;x ðxÞB I ðx
T
xI Þ þ b ðxÞB I;x ðx xI Þ
(29)
and the second derivative of the shape function can be calculated again by taking derivative of Eq. (29): cI;xx ðxÞ ¼ bT;xx ðxÞB I ðx xI Þ þ 2bT;x ðxÞB I;x ðx xI Þ þ bT ðxÞBI;xx ðx xI Þ.
ð30Þ
In this work, the cubic spline function is chosen as the weight function. Thus for this two-dimensional problem, the kernel function is expressed by Fa ðx xI Þ ¼ fx ðx xI Þ fy ðy yI Þ,
(31)
where weight function f is given by 9 82 2 3 for 0pjzI jp 12 > > = < 3 4zI þ 4zI fx ðzI Þ ¼ 43 4zI þ 4z2I 43 z3I for 12 ojzI jp1 , > > ; : 0 otherwise
also included b ¼ b¯
on l u ,
(36)
where dw , (37) dx and b¯ is the prescribed rotation on the boundary. The variational form due to the rotation is expressed as Z a¯ ¯ T ðb bÞ ¯ dl. Gb¯ ¼ ðb bÞ (38) 2 lu b¼
Therefore, the variational form due to the boundary conditions can be expressed as GB ¼ Gu¯ þ Gb¯ .
(39)
2.4. Instability regions via the kp-Ritz method and Bolotin’s method
ðx xI Þ , (32) a where the dilatation parameter, a, is the size of the support. At a node, the size of the domain of influence is calculated by
zI ¼
aI ¼ amax d I ,
1159
(33)
where amax is a scaling factor which ranges from 2.0 to 4.0. The distance d I is determined by searching for enough nodes to avoid singularity of the matrix M.
The total energy functional for this problem becomes G ¼ Gt þ GB .
(40)
Substituting the displacement functions of Eq. (15) into the total energy functional of Eq. (40) and applying the Rayleigh–Ritz minimization procedure with qG ¼ 0; qD
D ¼ uI ; vI ; wI ; I ¼ 1; 2; . . . ; NP
(41)
a system of Mathieu–Hill equations are obtained 2.3. Enforcement of boundary conditions—a penalty approach
~ q þ ðK ~ cos PtQÞq ¼ 0, M€
(42)
where It is necessary to deal with different boundary conditions to solve the problem. There are several approaches to enforce essential boundary conditions for meshless methods, such as Lagrange multipliers approach, penalty method approach, modified variational principles, etc. In the present work, the penalty method, see Reddy [26,27], is utilized to implement essential boundary conditions. The penalty formulation is developed as follows: Simply supported boundary conditions: For the domain bounded by l u , the displacement boundary condition is u ¼ u¯
on l u
(34)
in which u¯ is the prescribed displacement on the displacement boundary l u . The variational form is given by Z a¯ Gu¯ ¼ ðu u¯ ÞT ðu u¯ Þ dl, (35) 2 lu where a¯ is the penalty parameter taken as 103 E 11 , with E 11 being the elastic modulus of the shell in the principal coordinate direction. Clamped boundary conditions: In the clamped case, for the domain bounded by l u , besides the boundary condition described by Eq. (34), the rotation boundary condition is
~ ¼ K1 KKT ; K KIJ ¼ cI ðxJ ÞI;
~ ¼ K1 MK ¯ T , M I is the identity matrix.
A
B1
B2
K¼K þK þK þK Z L Z y0 KIJ ¼ R BI T SBJ dy dx, 0
u
KA IJ ¼ R
Z
L 0
ð44Þ ð45Þ
0
Z y0 ¯ M ¼ rhR MTI MJ dy dx 0 0 Z T KBIJ1 ¼ a¯ B1BI B1BJ dl, Zu T B2 KIJ ¼ a¯ B2BI B2BJ dl, Z
ð43Þ
L
Z
y0 0
T
A BA I N 0 BJ dy dx,
Z L Z y0 T A QA ¼ R BA IJ I N s BJ dy dx, 0 0 2 3 cI 0 0 6 7 0 cI 0 7, MI ¼ 6 4 5 0 0 cI
ð46Þ ð47Þ ð48Þ ð49Þ ð50Þ
ð51Þ
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1160
3
2
qcI 6 qx 6 6 6 6 0 6 6 6 1 qc 6 I 6 6 R qy 6 BI ¼ 6 6 6 0 6 6 6 6 6 0 6 6 6 4 0 2
cI 6 0 B1BI ¼ 6 4 0
0 1 qcI R qy qcI qx 0 1 qcI R2 qy 2 qcI R qx 0 cI 0
0
3. Numerical results and discussions
0
3
7 0 7; 5 cI
7 7 7 7 cI 7 7 R 7 7 7 0 7 7 7 , 2 q cI 7 7 2 7 qx 7 7 1 q2 c I 7 7 2 2 7 R qy 7 7 7 2 q2 c I 5 R qxqy 2 cI;x 6 0 B2BI ¼ 6 4 0
ð52Þ
N buc ¼
0 cI;x 0
0
3
7 0 7. ð53Þ 5 cI;x
Eq. (42) is in the form of a second order differential equation with periodic coefficients of the Mathieu–Hill type. The regions of instability are separated by the periodic solutions having periods T and 2T with T ¼ 2p=P. The solutions with period 2T are of greater practical importance as the widths of these unstable regions are usually larger than those associated with solutions having period of T. Using Bolotin’s approach, as a first approximation, the periodic solutions with period 2T can be sought in the form Pt Pt q¯ ¼ f¯ sin þ g¯ cos , 2 2
(54)
where f¯ and g¯ are arbitrary vectors. Substituting Eq. (54) into Eq. (42) and equating the coefficients of sin Pt=2 and cos Pt=2 terms, a set of linear homogeneous algebraic equations in terms of f¯ and g¯ can be obtained. The condition of non-trivial solutions are given by 20
det 4@ ¼ 0.
~ IJ þ K ~ IJ 1 QIJ 14 P2 M 2
0
0
~ IJ þ K ~ IJ þ 1 QIJ 14 P2 M 2
13
Eh2 ½3ð1 n2 Þ1=2 R
,
(57)
where E is the elastic modulus and n is the Poisson’s ratio of the isotropic cylindrical shell. For laminated circular cylindrical panels, the critical buckling load N cr is approximated as N cr ¼
E 22 h2 ½3ð1 n12 n21 Þ1=2 R
.
(58)
In order to evaluate the validity and accuracy of the present formulation, the points of origin of first four unstable regions for a panel having simply supported boundary condition under tensile and compressive loadings are compared with those given by Ng et al. [15]. The panel has parameters of b=R ¼ 0:5, L=R ¼ 2 and n ¼ 0:3 The non-dimensional frequency parameter P¯ is defined as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ¯ ¼ RP rð1 v Þ. (59) P E Each unstable regions is bounded by two curves originating from a common point from the P¯ axis with N s ¼ 0. The numerical results are tabulated in Table 1. The thickness ratio of h=R ranges from 0.03 to 0.05. The first four modes are ðm ¼ 1; n ¼ 1Þ, ðm ¼ 2; n ¼ 1Þ, ðm ¼ 3; n ¼ 1Þ, and ðm ¼ 4; n ¼ 1Þ. It is evident that the present results agree well with those reported by Ng et al. [15] for the panel under the tensile loading. The present solutions are a little higher for the panel under the compressive loading. The maximum difference is below 4%.
A5
3.1. Isotropic cylindrical panels ð55Þ
Instead of solving the above nonlinear geometric equations for P, the above expression can be rearranged to the standard form of a generalized eigenvalue problem 20 1 0 13 1 ~ ~ IJ 1 QIJ K 0 0 MIJ 2 4 A P2 @ A5 det4@ 1 ~ ~ IJ þ 1 QIJ M 0 K 0 IJ 2 4 ¼ 0.
To assess the stability and accuracy of the present methodology, numerical comparisons and convergence studies are performed. For periodic compressive loads, it is obvious that the compressive axial loads cannot exceed the critical buckling load of the cylindrical shell panel. In this paper, for isotropic cylindrical shell panel of intermediate length, the buckling load is given in Timoshenko and Gere [28] as
ð56Þ
The generalized eigenvalues P2 of the above generalized eigenvalue problem define the boundaries between the stable and unstable regions.
The instability regions of isotropic cylindrical panels having different boundary conditions thickness ratios are investigated here. The panel parameters are chosen to be y0 ¼ 60 , L=R ¼ 2 and n ¼ 0:3. The thickness ratio, h=R, 1 1 varies from 100 to 300 . The static axial force N 0 ¼ 0:5N buc is selected. Fig. 2 shows the first four unstable regions of a 1 panel having h=R ¼ 100 under simply supported boundary conditions. These first four modes, respectively, correspond to ðm ¼ 2; n ¼ 1Þ, ðm ¼ 1; n ¼ 1Þ, ðm ¼ 3; n ¼ 1Þ, and ðm ¼ 2; n ¼ 2Þ. It is observed that the fourth mode ðm ¼ 2; n ¼ 2Þ, has the largest instability region, followed in descending order by modes ðm ¼ 2; n ¼ 1Þ, ðm ¼ 1; n ¼ 1Þ, and ðm ¼ 3; n ¼ 1Þ, which has the smallest
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
0.50
Table 1 Comparison of unstable region for a panel having simply supported boundary condition under tensile and compressive loading ðb=R ¼ 0:5; L=R ¼ 2:0; N 0 ¼ 0:5N cr Þ Mode
0.45
Point of origin Tensile loading
0.03
0.04
0.05
Compressive loading
Ng et al. [15]
Present
Ng et al. [15]
Present
(1,1) (2,1) (3,1) (4,1)
0.7651 1.0722 1.5076 1.9976
0.7512 1.0661 1.4929 1.9973
0.6533 0.7182 0.9172 1.1958
0.6352 0.6956 0.8564 1.1371
(1,1) (2,1) (3,1) (4,1)
0.9977 1.3375 1.8313 2.4111
0.9756 1.3229 1.8140 2.4148
0.8854 0.9716 1.2024 1.5495
0.8839 0.9520 1.1409 1.4969
(1,1) (2,1) (3,1) (4,1)
1.2308 1.6063 2.1621 2.8345
1.2008 1.5833 2.1408 2.8405
1.1183 1.2347 1.5135 1.9419
1.0836 1.1921 1.4489 1.8929
0.40
mode 1 P
h=R
1161
0.35
mode 2 mode 3 mode 4
0.30
0.25
0.20 0.65
0.0
0.1
0.60
0.2
0.3 Ns / No
0.4
0.5
Fig. 3. Unstable regions for the first four modes of cylindrical panel with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 200Þ.
0.55
P
0.50
0.34
0.45
0.32 0.40 mode 1
mode 3
mode 2
mode 4
0.30
0.35
0.28 0.30
mode 1
0.26 P
mode 2
0.25 0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 2. Unstable regions for the first four modes of cylindrical panel with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
instability region. It is also noted that the region size associated with the fourth mode ðm ¼ 2; n ¼ 2Þ is approximately twice as wide as the principal region of the fundamental mode ðm ¼ 2; n ¼ 1Þ. Fig. 3 presents the first four instability regions of a simply supported shell panel 1 with thickness ratio h=R ¼ 200 . The first four modes are, respectively, modes ðm ¼ 1; n ¼ 2Þ, ðm ¼ 1; n ¼ 3Þ,
mode 3
0.24
mode 4 0.22 0.20 0.18 0.16 0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 4. Unstable regions for the first four modes of cylindrical panel with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 300Þ.
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1162
0.40
0.14 0.13
0.35
0.12 0.30
0.11 0.10 P
P
mode 1 0.25
mode 2
mode 1 mode 2 mode 3 mode 4
0.09
mode 3 mode 4 0.20
0.08 0.07
0.15
0.06 0.10
0.05 0.0
0.0
0.1
0.2 Ns / No
0.3
0.4
0.1
0.5
0.2 Ns / No
0.3
0.4
0.5
Fig. 7. Unstable regions for the first four modes of cylindrical panel with CFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 300Þ.
Fig. 5. Unstable regions for the first four modes of cylindrical panel with CFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
0.26 0.18
0.24 0.16
mode 1
mode 1 mode 2 mode 3 mode 4
0.22
mode 2 mode 3
0.14
0.20 P
P
mode 4 0.12
0.18 0.10
0.16 0.08
0.14 0.06 0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 6. Unstable regions for the first four modes of cylindrical panel with CFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 200Þ.
0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 8. Unstable regions for the first four modes of cylindrical panel with FFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1163
0.71
0.16 mode 1 mode 2 mode 3 mode 4
0.15
o
o
o
0 /90 /90 /0
o
o
o
o
o
0 /90 /0 /90
0.70
ϖ
0.14 0.69
P
0.13
0.12
0.68
0.11 0.67 0.0
0.1
0.2
0.3
0.4
0.5
Ns / No
0.10
0.09 0.0
0.1
0.2
0.3
0.4
0.5
Fig. 11. Unstable regions for (0 =90 =0 =90 and 0 =90 =90 =0 ) cylindrical panels with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
Ns / No 0.89
Fig. 9. Unstable regions for the first four modes of cylindrical panel with FFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 200Þ.
0.88 0.87
0.12
45o/-45o 45o/-45o/45o/-45o 45o/-45o/45o/-45o/45o/-45o
0.86 ϖ
0.13 mode 1 mode 2 mode 3 mode 4
0.85 0.84 0.83
0.11
P
0.82 0.0
0.1
0.2
0.3 Ns / No
0.4
0.5
0.10 Fig. 12. Unstable regions for simply supportive cylindrical panels with different layers ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
0.09
0.08
0.0
0.1
0.2
0.3
0.4
0.5
Ns / No Fig. 10. Unstable regions for the first four modes of cylindrical panel with FFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 300Þ.
ðm ¼ 1; n ¼ 1Þ, and ðm ¼ 1; n ¼ 4Þ. It is clear that the first four modes in Fig. 3 are of different order as with those in Fig. 2 due to the change in panel thickness. It is also observed that the widths of the instability regions of panels are narrower as compared with the corresponding h=R ¼ 1 100 cases, and increase in the sequence of modes ðm ¼ 1; n ¼ 2Þ, ðm ¼ 1; n ¼ 3Þ, ðm ¼ 1; n ¼ 1Þ, andðm ¼ 1; n ¼ 4Þ. It is further noted that the sizes of the instability regions of the two lowest modes, ðm ¼ 1; n ¼ 1Þ and ðm ¼ 1; n ¼ 4Þ, have no significant difference. For a panel 1 with thickness ratio h=R ¼ 300 , corresponding results are
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1164
frequency parameter o ¯ is defined as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rð1 n12 n21 Þ . o ¯ ¼ RP E 22
0.72 α = 15o α = 30o
(60)
o
α = 45
0.68
α = 60o o
ϖ
α = 75 0.64
0.60
0.56 0.0
0.1
0.2
0.3
0.4
0.5
Ns / No Fig. 13. Unstable regions for ða= a=a= aÞ cylindrical panels with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ. 1 shown in Fig. 4, and similar trends to that of the h=R ¼ 100 case are observed. However, the first four modes are now ðm ¼ 1; n ¼ 2Þ, ðm ¼ 1; n ¼ 3Þ, ðm ¼ 1; n ¼ 4Þ, and ðm ¼ 2; n ¼ 3Þ, respectively. Figs. 5–7 show the first four principal instability regions of the CFFF cylindrical panels with respective thickness 1 1 1 ratios of h=R ¼ 100 , h=R ¼ 200 , and h=R ¼ 300 . For the panel 1 having thickness ratio h=R ¼ 100 , the instability region of mode 4 is the widest and much larger than that of modes 1–3. The instability region for mode 2 is larger than that of mode 3 but less than that of mode 1. Further, it is noticed that the instability regions of modes 1 and 2 are overlapping as N s =N 0 exceeds 0.33. The first four unstable 1 regions of a panel with h=R ¼ 200 show different trends. From Fig. 6, it is observed that the width of the unstable regions reduces from mode 1 to 4, and some portions of the unstable regions of modes 1 and 2 are overlapping. This trend is also observed from Fig. 7 for the panel having 1 thickness ratio h=R ¼ 300 . The respective unstable regions of completely free panels 1 1 1 with thickness ratios h=R ¼ 100 , h=R ¼ 200 , and h=R ¼ 300 are depicted in Figs. 8–10. It is obvious that all three figures show similar general trends, and the widths of the instability regions for all four modes in each case have no significant differences. It is also found that certain portions of the unstable regions of all four modes mutually overlap 1 for the panels with thickness ratios h=R ¼ 200 and 1 h=R ¼ 300.
3.2. Laminated composite cylindrical panels The material properties of the laminated cylindrical panels in the present study are E 1 =E 2 ¼ 40, G 12 =E 2 ¼ 0:5, h=R ¼ 1=100, and n ¼ 0:25. The static axial compressive force is chosen to be N 0 ¼ 0:3N cr . The non-dimensional
For the simply supported panels with the symmetric and anti-symmetric cross-ply lamination schemes, 0 =90 =0 = 90 and 0 =90 =90 =0 , the unstable regions are plotted in Fig. 11. The modes for the two lamination schemes correspond to ðm ¼ 1; n ¼ 1Þ. It is observed that the two regions have comparatively same width, and the instability region for panel with 0 =90 =0 =90 scheme is at higher frequencies than that of panel with 0 =90 =90 =0 scheme. The effects of number of layers on the instability regions of laminated panels are demonstrated in Fig. 12. The twolayer, four-layer and six-layer panels, which have lamination schemes of 45 = 45 , 45 = 45 =45 = 45 , and 45 = 45 =45 = 45 =45 = 45 , respectively, are investigated. The corresponding mode in this case is ðm ¼ 1; n ¼ 1Þ. It is noted that all instability regions have higher frequencies and become narrower as the number of layers increases. It is also seen that the instability regions for four-layer and six-layer laminates begin to overlap as N s =N 0 exceeds 0.31. It is expected that the overlapping area of instability regions will become larger when the number of layers increases further. Fig. 13 illustrates the influence of lamination angles on the instability regions of laminates. The lamination scheme a= a=a= a is selected, five lamination angles, 15 , 30 , 45 , 60 , and 75 , are considered in this case. It is observed that, except mode ðm ¼ 1; n ¼ 1Þ corresponded by panel with lamination angle 60 , the mode ðm ¼ 1; n ¼ 2Þ, relates to the panels with other four lamination angles. It is also noted that the instability region of panel with angle 45 has the highest frequencies, followed by those of panels with lamination angles of 15 , 60 , 30 , and 75 . The instability region widens in the same lamination angle sequence.
4. Conclusions The dynamic stability of thin circular cylindrical panels under static and periodic axial forces has been investigated using the mesh-free kp-Ritz method. The mesh-free kernel particle estimate was successfully employed to approximate the 2D transverse displacement field. The Ritz minimization procedure was employed to obtain a system of Mathieu–Hill equations, and the principal instability regions were calculated via Bolotin’s first approximation. Numerical comparisons were made with those in literature, and the present methodology validated. The instability regions of composite laminated cylindrical panels having different lamination schemes and layers were investigated. It was found that the both the lamination scheme and number of layer influenced the positions and sizes of the instability regions.
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
Acknowledgements The work described in this paper was fully supported by a grant from the Research Grants Councils of the Hong Kong Special Administrative Region, China (9041011 CityU 1160/05E). References [1] Bolotin VV. The dynamic stability of elastic systems. Holden-Day, San Francisco; 1964. [2] Vijayaraghavan A, Evan-Iwanowski RM. Parametric instability of circular cylindrical shells. TASME, Journal of Applied Mechanics 1967;31:985–90. [3] Yao JC. Nonlinear elastic buckling and parametric excitation of a cylinder under axial loads. TASME, Journal of Applied Mechanics 1965;29:109–15. [4] Birman V. Dynamic stability of unsymmetrically laminated rectangular plates. Mechanics Research Communications 1985;12:81–6. [5] Srinivasan RS, Chellapandi P. Dynamic stability of rectangular laminated composite plates. Computers and Structures 1986;24:233–8. [6] Moorthy J, Reddy JN. Parametric instability of laminated composite plates with transverse shear deformation. International Journal of Solids and Structures 1990;26:801–11. [7] Bert CW, Birman V. Dynamic stability of shear deformable antisymmetric angle-ply plates. International Journal of Solids and Structures 1987;23:1053–61. [8] Bert CW, Birman V. Parametric instability of thick, orthotropic, circular cylindrical shells. Acta Mechanica 1988;71:61–76. [9] Koval LR. Effects of longitudinal resonance on the parametric stability of an axially excited cylindrical shell. Journal of Acoustical Society of America 1974;55:91–7. [10] Argento A, Scott RA. Dynamic instability of laminated anisotropic circular cylindrical shells, part I: theoretical development. Journal of Sound and Vibration 1993;162:311–22. [11] Argento A, Scott RA. Dynamic instability of laminated anisotropic circular cylindrical shells, part II: numerical results. Journal of Sound and Vibration 1993;162:323–32. [12] Ng TY, Lam KY, Reddy JN. Dynamic stability of cross-ply laminated composite cylindrical shells, International. Journal of Mechanical Sciences 1998;40:805–23. [13] Lam KY, Ng TY. Dynamic stability of cylindrical shells subjected to conservative periodic axial loads using different shell theories. Journal of Sound and Vibration 1997;207:497–520.
1165
[14] Ganapathi M, Varadan TK, Balamurugan V. Dynamic instability of laminated composite curved panels using finite element method. Computers and Structures 1994;53:335–42. [15] Ng TY, Lam KY, Reddy JN. Dynamic stability of cylindrical panels with transverse shear effects. International Journal of Solids and Structures 1999;36:3483–96. [16] Kitipornchai S, Xiang Y, Wang CM, Liew KM. Buckling of thick skew plates. International Journal for Numerical Method in Engineering 1993;36:1299–310. [17] Liew KM, Xiang Y, Wang CM, Kitipornchai S. Flexural vibration of shear deformable circular and annular plates on ring supports. Computer Methods in Applied Mechanics and Engineering 1993;110:301–15. [18] Liew KM, Xiang Y, Kitipornchai S, Meek JL. Formulation of Mindlin-Engesser model for stiffened plate vibration. Computer Methods in Applied Mechanics and Engineering 1995;120: 339–53. [19] Cheung YK, Zhou D. Free vibrations of rectangular unsymmetrically laminated composite plates with internal line supports. Computers and Structures 2001;79:1923–32. [20] Cheung YK, Zhou D. Three-dimensional vibration analysis of cantilevered and completely free isosceles triangular plates. International Journal of Solids Structures 2002;39:673–87. [21] Zeng H, Bert CW. Free vibration analysis of discretely stiffened skew plates. International Journal of Stability and Dynamics 2001;1: 125–44. [22] Su GH, Xiang Y. A non-discrete approach for analysis of plates with multiple subdomains. Engineering Structures 2002;24: 563–75. [23] Lim CW, Liew KM. A pb-2 Ritz formulation for flexural vibration of shallow cylindrical shells of rectangular planform. Journal of Sound and Vibration 1994;173:343–75. [24] Liew KM, Feng ZC. Three-dimensional free vibration analysis of perforated super elliptical plates via the p-Ritz method. International Journal of Mechanical Sciences 2001;43:2613–30. [25] Chen JS, Pan C, Wu CT. Large deformation analysis of rubber based on a reproducing kernel particle method. Computational Mechanics 1997;19:211–27. [26] Reddy JN. An introduction to the finite element method. New York: McGraw-Hill; 1993. [27] Reddy JN. Applied functional analysis and variational methods in engineering. New York: McGraw-Hill; 1986; reprinted by Krieger, Melbourne, Florida, 1991. [28] Timoshenko SP, Gere JM. Theory of elastic stability. New York: McGraw-Hill; 1961.
International Journal of Mechanical Sciences 49 (2007) 1156–1165 www.elsevier.com/locate/ijmecsci
Dynamic stability analysis of composite laminated cylindrical panels via the mesh-free kp-Ritz method K.M. Liewa,, Y.Y. Leea, T.Y. Ngb, X. Zhaob a
b
Department of Building and Construction, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong School of Mechanical and Production Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore Received 14 June 2006; received in revised form 19 February 2007; accepted 19 February 2007 Available online 24 February 2007
Abstract In this paper, the dynamic stability analysis of thin, laminated cylindrical panels under static and periodic axial forces is presented by using the mesh-free kp-Ritz method. The mesh-free kernel particle estimate is employed to approximate the 2D transverse displacement field. A system of Mathieu–Hill equations is obtained through the application of the Ritz minimization procedure to the energy expressions. The principal instability regions are then analyzed via Bolotin’s first approximation. Effects of lamination schemes such as number of layers, ply-angle and boundary conditions on the instability regions are examined in detail. r 2007 Elsevier Ltd. All rights reserved. Keywords: Dynamic stability; Parametric resonance; Cylindrical panel; Composite laminate; Mesh-free kp-Ritz method; Bolotin’s first approximation
1. Introduction The dynamic stability or phenomenon of parametric resonance of a structure under periodic loads has attracted an enormous research interest. The earlier work include the dynamic stability analysis of elastic system by Bolotin [1], the parametric instability of circular cylindrical shells by Vijayaraghavan and Evan-Iwanowski [2], and nonlinear elastic buckling and parametric excitation of a cylinder under axial loads by Yao [3]. Birman [4] and Srinivasan and Chellapandi [5] studied the dynamic stability of laminated rectangular plates. The parametric instability of plates with transverse shear deformation was investigated by Moorthy and Reddy [6] and Bert and Birman [7]. Bert and Birman [8] extended Yao’s [3] approach to the parametric instability of thick orthotropic shells using higher-order theory. Koval [9] analyzed the effects of longitudinal resonance on the parametric stability of an axially excited cylindrical shell. Argento and Scott [10,11] employed a perturbation technique to examine the dynamic stability of layered anisotropic circular cylindrical shells Corresponding author. Tel.: +852 3442 65816; fax: +852 2788 7612.
E-mail address: [email protected] (K.M. Liew). 0020-7403/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2007.02.005
under axial loading. More recently, the dynamic stability of cross-ply laminated composite cylindrical shells under combined static and periodical axial forces was investigated by Ng et al. [12] and Lam and Ng [13]. The dynamic instability of laminated composite curved panels was presented by Ganapathi et al. [14] using finite element method. Ng et al. [15] investigated the dynamic stability of a simply supported cylindrical panel using generalizes Donnell’s shell theory. The Ritz or Rayleigh–Ritz method is a proven approximate technique in computational mechanics, and notable works include those of Kitipornchai et al. [16], Liew et al. [17,18], Cheung and Zhou [19,20], Zeng and Bert [21], Su and Xiang [22], Lim and Liew [23], and Liew and Feng [24]. In this paper, for the dynamic stability analysis of composite laminated cylindrical shell panels under combined static and periodic axial loads, an energy formulation is first described, and a system of Mathieu–Hill equations is obtained via the Ritz minimization procedure. The parametric resonance responses are then analyzed using Bolotin’s [1] method. In conventional Ritz method, it is difficult to choose appropriate functions to satisfy some complicate boundary conditions, and the eigen equations need to be reevaluated for different boundary condition
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
cases. In the present kp-Ritz method, a standard weight function is employed to express the interior field, and the boundary conditions are enforced by penalty method. Therefore, the disadvantages in conventional Ritz method are avoided. For the assumed functions in the present work, the mesh-free kernel particle approximate is utilized to describe the 2D transverse displacement field. The commonly used trigonometric or hyperbolic functions in conjunction with the Ritz procedure are not suitable for general laminates as erroneously vanish. The present formulation ensures that the bending–extension coupling stiffnesses in general laminates are retained. The present method is verified by comparing with those results in open literature. The effects of boundary conditions and lamination schemes on the instability regions are also investigated. 2. Theoretical formulation
1157
load is given by N a ¼ N 0 þ N s cos Pt,
(1)
where P is the frequency of excitation in radians per unit time. The kinetic energy for the cylindrical panels can be expressed as Z L Z y0 1 ½u_ 2 þ v_2 þ w_ 2 R dy dx, (2) T ¼ rh 2 0 0 where the three terms are due to contributions from the linear velocities in the x, y, and z directions, respectively. The strain energy U a due to the axial loading can be written as " 2 2 # Z Z 1 L y0 qu 2 qv qw Ua ¼ Na þ þ R dy dx. 2 0 0 qx qx qx (3)
2.1. Energy formulation for cylindrical panels The cylindrical panel considered is shown in Fig. 1a, where a coordinate system ðx; y; zÞ is fixed on the middle surface of the panel. This panel is considered to be thin and of length L, radius R and thickness h, bounded along its edges by the lines x ¼ 0, x ¼ L, y ¼ 0 and y ¼ y0 . The displacements of the panel in the x, y, and z directions are denoted by u, v and w, respectively. The pulsating axial
s
L
w
h
v θο
b
u
R θ
T ¼ fe1 e2 g k1 k2 2tg,
(5)
where the middle surface strains, e1 , e2 and g, and the middle surface curvatures, k1 , k2 and t, are defined according to Love’s thin shell theory qu 1 qv qv 1 qu ; e2 ¼ þw ; g¼ þ e1 ¼ qx R qy qx R qy
and [S] is given by 2 A11 A12 A16 6 6 A12 A22 A26 6 6 A16 A26 A66 6 ½S ¼ 6 6 B11 B12 B16 6 6 B12 B22 B26 4 B16 B26 B66
v
h2
h4 h5
h
B11
B12
B16
B12 B16
B22 B26
B26 B66
D11
D12
D12 D16
D22 D26
ð6Þ
3
7 7 7 7 7 7, D16 7 7 D26 7 5 D66
(7)
where the extensional ðAij Þ, coupling ðBij Þ and bending ðDij Þ stiffnesses, are defined as Z
Fig. 1. (a) Coordinate system of the cylindrical panel; (b) cross-sectional view of the laminated cylindrical panel.
(4)
where T and [S] are the strain vector and stiffness matrix, respectively, and T can be defined as
q2 w 1 q2 w qv k1 ¼ 2 ; k2 ¼ 2 , R q2 y qy q x 1 q2 w qv t¼ R qxqy qx
w
h1
The strain energy of the shell is given by Z Z 1 L y0 T U ¼ ½SR dy dx, 2 0 0
h=2
ðAij ; Bij ; Dij Þ ¼ h=2
Qij ð1; z; z2 Þ dz.
(8)
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1158
For a panel composed of different layers of materials, the stiffnesses are expressed as Nl X ¯ kij ðhk hkþ1 Þ, Q Aij ¼ k¼1
Bij ¼
l 1X ¯ k ðh2 h2kþ1 Þ, Q 2 k¼1 ij k
Dij ¼
Nl 1X ¯ k ðh3 h3 Þ, Q kþ1 3 k¼1 ij k
N
cI ðxÞ ¼ Cðx; x xI ÞFa ðx xI Þ,
ð9Þ
Cðx; x xI Þ ¼ H T ðx xI ÞbðxÞ,
cos2 a sin2 a
where a is the angular orientation of the fibres and [Q] is the reduced stiffness matrix defined as 2 3 Q11 Q12 0 6 0 7 (12) ½Q ¼ 4 Q12 Q22 5. 0 0 Q66 The material constants in the reduced stiffness matrix [Q] are given as E 11 n12 E 22 Q11 ¼ ; Q12 ¼ , 1 n12 n21 1 n12 n21 E 22 Q22 ¼ ; Q66 ¼ G 12 , ð13Þ 1 n12 n21 where E 11 and E 22 are the elastic moduli in the principle material coordinates, G 12 is the shear modulus and n12 and n21 are Poisson’s ratios. The total energy functional of the panel is thus Gt ¼ T U a U . (14)
bT ðxÞ ¼ ½b0 ðx; yÞ; b1 ðx; yÞ; b2 ðx; yÞ; b3 ðx; yÞ; b4 ðx; yÞ; b5 ðx; yÞ, H T ðx xI Þ ¼ ½1; x xI ; y yI ; ðx xI Þðy yI Þ, ðx xI Þ2 ; ðy yI Þ2 .
2.2. Two-dimensional kernel particle shape functions The discrete displacement approximations for shell panels can be expressed as
cI ðxÞ ¼ bT ðxÞHðx xI ÞFa ðx xI Þ. cI ðxÞ ¼ bT ðxÞB I ðx xI Þ,
vðx; yÞ ¼
where bðxÞ ¼ M 1 ðxÞHð0Þ,
ð21Þ
BI ðx xI Þ ¼ Hðx xI ÞFa ðx xI Þ
ð22Þ
and the moment matrix M is a function of x while Hð0Þ is a constant vector. The expressions of M and Hð0Þ are given by MðxÞ ¼
NP X
Hðx xI ÞH T ðx xÞFa ðx xI Þ,
NP X I¼1
cI ðx; yÞwI eiot ,
ð23Þ
I¼1
Hð0Þ ¼ ½1; 0; 0; 0; 0; 0T .
ð24Þ
The shape function can therefore be expressed as cI ðxÞ ¼ H T ð0ÞM 1 ðxÞHðx xI ÞFa ðx xI Þ.
(25)
For this thin shell problem, the first and second derivatives of shape functions need to be determined. Eq. (21) can be rewritten as (26)
The vector bðxÞ can be determined using LU decomposition of the matrix MðxÞ followed by backward substitution. The derivatives of bðxÞ can be obtained similarly. Taking the derivative of Eq. (26), we have
MðxÞb;x ðxÞ ¼ H ;x ð0Þ M ;x ðxÞbðxÞ.
cI ðx; yÞvI eiot ,
I¼1
wðx; yÞ ¼
(20)
(27)
which can be rearranged as
I¼1 NP X
(19)
Eq. (16) can be rewritten as
M ;x ðxÞbðxÞ þ MðxÞb;x ðxÞ ¼ H ;x ð0Þ
cI ðx; yÞuI eiot ,
ð18Þ
H is a vector of quadratic basis and bi ðx; yÞ’s are functions of x and y which are to be determined. Thus, the shape function can be written as
MðxÞbðxÞ ¼ Hð0Þ.
NP X
(17)
where
where [T] is the transformation matrix for the principle material coordinates and the panel’s coordinates, and is defined as 2 3 sin2 a 2 cos a sin a cos2 a 6 7 ½T ¼ 4 sin2 a (11) cos2 a 2 cos a sin a 5,
uðx; yÞ ¼
(16)
where Cðx; x xI Þ is the correction function and Fa ðx xI Þ is the kernel function. The correction function Cðx; x xI Þ is described as
where hk and hkþ1 denote the distances from the panel reference surface to the outer and inner surfaces of the kth layer, respectively, as shown in Fig. 1b. N l denotes the total ¯ kij is the number of layers in the laminated panel and Q transformed reduced stiffness matrix for the kth layer defined as ¯ ¼ ½T½Q½TT , ½Q (10)
cos a sin a cos a sin a
where NP is the total number of particles and cI ðx; yÞ’s are the shape functions of displacements u, v, and w. The 2D shape function is expressed as (see [25])
ð15Þ
(28)
It can be seen that the first derivative of bðxÞ can be formulated using the LU decomposition procedure again. The second derivative of bðxÞ can be determined by taking derivative of Eq. (28) and using the same LU
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
decomposition procedure. The first derivative of the shape function can be obtained by taking the derivative of Eq. (20): cI;x ðxÞ ¼
bT;x ðxÞB I ðx
T
xI Þ þ b ðxÞB I;x ðx xI Þ
(29)
and the second derivative of the shape function can be calculated again by taking derivative of Eq. (29): cI;xx ðxÞ ¼ bT;xx ðxÞB I ðx xI Þ þ 2bT;x ðxÞB I;x ðx xI Þ þ bT ðxÞBI;xx ðx xI Þ.
ð30Þ
In this work, the cubic spline function is chosen as the weight function. Thus for this two-dimensional problem, the kernel function is expressed by Fa ðx xI Þ ¼ fx ðx xI Þ fy ðy yI Þ,
(31)
where weight function f is given by 9 82 2 3 for 0pjzI jp 12 > > = < 3 4zI þ 4zI fx ðzI Þ ¼ 43 4zI þ 4z2I 43 z3I for 12 ojzI jp1 , > > ; : 0 otherwise
also included b ¼ b¯
on l u ,
(36)
where dw , (37) dx and b¯ is the prescribed rotation on the boundary. The variational form due to the rotation is expressed as Z a¯ ¯ T ðb bÞ ¯ dl. Gb¯ ¼ ðb bÞ (38) 2 lu b¼
Therefore, the variational form due to the boundary conditions can be expressed as GB ¼ Gu¯ þ Gb¯ .
(39)
2.4. Instability regions via the kp-Ritz method and Bolotin’s method
ðx xI Þ , (32) a where the dilatation parameter, a, is the size of the support. At a node, the size of the domain of influence is calculated by
zI ¼
aI ¼ amax d I ,
1159
(33)
where amax is a scaling factor which ranges from 2.0 to 4.0. The distance d I is determined by searching for enough nodes to avoid singularity of the matrix M.
The total energy functional for this problem becomes G ¼ Gt þ GB .
(40)
Substituting the displacement functions of Eq. (15) into the total energy functional of Eq. (40) and applying the Rayleigh–Ritz minimization procedure with qG ¼ 0; qD
D ¼ uI ; vI ; wI ; I ¼ 1; 2; . . . ; NP
(41)
a system of Mathieu–Hill equations are obtained 2.3. Enforcement of boundary conditions—a penalty approach
~ q þ ðK ~ cos PtQÞq ¼ 0, M€
(42)
where It is necessary to deal with different boundary conditions to solve the problem. There are several approaches to enforce essential boundary conditions for meshless methods, such as Lagrange multipliers approach, penalty method approach, modified variational principles, etc. In the present work, the penalty method, see Reddy [26,27], is utilized to implement essential boundary conditions. The penalty formulation is developed as follows: Simply supported boundary conditions: For the domain bounded by l u , the displacement boundary condition is u ¼ u¯
on l u
(34)
in which u¯ is the prescribed displacement on the displacement boundary l u . The variational form is given by Z a¯ Gu¯ ¼ ðu u¯ ÞT ðu u¯ Þ dl, (35) 2 lu where a¯ is the penalty parameter taken as 103 E 11 , with E 11 being the elastic modulus of the shell in the principal coordinate direction. Clamped boundary conditions: In the clamped case, for the domain bounded by l u , besides the boundary condition described by Eq. (34), the rotation boundary condition is
~ ¼ K1 KKT ; K KIJ ¼ cI ðxJ ÞI;
~ ¼ K1 MK ¯ T , M I is the identity matrix.
A
B1
B2
K¼K þK þK þK Z L Z y0 KIJ ¼ R BI T SBJ dy dx, 0
u
KA IJ ¼ R
Z
L 0
ð44Þ ð45Þ
0
Z y0 ¯ M ¼ rhR MTI MJ dy dx 0 0 Z T KBIJ1 ¼ a¯ B1BI B1BJ dl, Zu T B2 KIJ ¼ a¯ B2BI B2BJ dl, Z
ð43Þ
L
Z
y0 0
T
A BA I N 0 BJ dy dx,
Z L Z y0 T A QA ¼ R BA IJ I N s BJ dy dx, 0 0 2 3 cI 0 0 6 7 0 cI 0 7, MI ¼ 6 4 5 0 0 cI
ð46Þ ð47Þ ð48Þ ð49Þ ð50Þ
ð51Þ
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1160
3
2
qcI 6 qx 6 6 6 6 0 6 6 6 1 qc 6 I 6 6 R qy 6 BI ¼ 6 6 6 0 6 6 6 6 6 0 6 6 6 4 0 2
cI 6 0 B1BI ¼ 6 4 0
0 1 qcI R qy qcI qx 0 1 qcI R2 qy 2 qcI R qx 0 cI 0
0
3. Numerical results and discussions
0
3
7 0 7; 5 cI
7 7 7 7 cI 7 7 R 7 7 7 0 7 7 7 , 2 q cI 7 7 2 7 qx 7 7 1 q2 c I 7 7 2 2 7 R qy 7 7 7 2 q2 c I 5 R qxqy 2 cI;x 6 0 B2BI ¼ 6 4 0
ð52Þ
N buc ¼
0 cI;x 0
0
3
7 0 7. ð53Þ 5 cI;x
Eq. (42) is in the form of a second order differential equation with periodic coefficients of the Mathieu–Hill type. The regions of instability are separated by the periodic solutions having periods T and 2T with T ¼ 2p=P. The solutions with period 2T are of greater practical importance as the widths of these unstable regions are usually larger than those associated with solutions having period of T. Using Bolotin’s approach, as a first approximation, the periodic solutions with period 2T can be sought in the form Pt Pt q¯ ¼ f¯ sin þ g¯ cos , 2 2
(54)
where f¯ and g¯ are arbitrary vectors. Substituting Eq. (54) into Eq. (42) and equating the coefficients of sin Pt=2 and cos Pt=2 terms, a set of linear homogeneous algebraic equations in terms of f¯ and g¯ can be obtained. The condition of non-trivial solutions are given by 20
det 4@ ¼ 0.
~ IJ þ K ~ IJ 1 QIJ 14 P2 M 2
0
0
~ IJ þ K ~ IJ þ 1 QIJ 14 P2 M 2
13
Eh2 ½3ð1 n2 Þ1=2 R
,
(57)
where E is the elastic modulus and n is the Poisson’s ratio of the isotropic cylindrical shell. For laminated circular cylindrical panels, the critical buckling load N cr is approximated as N cr ¼
E 22 h2 ½3ð1 n12 n21 Þ1=2 R
.
(58)
In order to evaluate the validity and accuracy of the present formulation, the points of origin of first four unstable regions for a panel having simply supported boundary condition under tensile and compressive loadings are compared with those given by Ng et al. [15]. The panel has parameters of b=R ¼ 0:5, L=R ¼ 2 and n ¼ 0:3 The non-dimensional frequency parameter P¯ is defined as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ¯ ¼ RP rð1 v Þ. (59) P E Each unstable regions is bounded by two curves originating from a common point from the P¯ axis with N s ¼ 0. The numerical results are tabulated in Table 1. The thickness ratio of h=R ranges from 0.03 to 0.05. The first four modes are ðm ¼ 1; n ¼ 1Þ, ðm ¼ 2; n ¼ 1Þ, ðm ¼ 3; n ¼ 1Þ, and ðm ¼ 4; n ¼ 1Þ. It is evident that the present results agree well with those reported by Ng et al. [15] for the panel under the tensile loading. The present solutions are a little higher for the panel under the compressive loading. The maximum difference is below 4%.
A5
3.1. Isotropic cylindrical panels ð55Þ
Instead of solving the above nonlinear geometric equations for P, the above expression can be rearranged to the standard form of a generalized eigenvalue problem 20 1 0 13 1 ~ ~ IJ 1 QIJ K 0 0 MIJ 2 4 A P2 @ A5 det4@ 1 ~ ~ IJ þ 1 QIJ M 0 K 0 IJ 2 4 ¼ 0.
To assess the stability and accuracy of the present methodology, numerical comparisons and convergence studies are performed. For periodic compressive loads, it is obvious that the compressive axial loads cannot exceed the critical buckling load of the cylindrical shell panel. In this paper, for isotropic cylindrical shell panel of intermediate length, the buckling load is given in Timoshenko and Gere [28] as
ð56Þ
The generalized eigenvalues P2 of the above generalized eigenvalue problem define the boundaries between the stable and unstable regions.
The instability regions of isotropic cylindrical panels having different boundary conditions thickness ratios are investigated here. The panel parameters are chosen to be y0 ¼ 60 , L=R ¼ 2 and n ¼ 0:3. The thickness ratio, h=R, 1 1 varies from 100 to 300 . The static axial force N 0 ¼ 0:5N buc is selected. Fig. 2 shows the first four unstable regions of a 1 panel having h=R ¼ 100 under simply supported boundary conditions. These first four modes, respectively, correspond to ðm ¼ 2; n ¼ 1Þ, ðm ¼ 1; n ¼ 1Þ, ðm ¼ 3; n ¼ 1Þ, and ðm ¼ 2; n ¼ 2Þ. It is observed that the fourth mode ðm ¼ 2; n ¼ 2Þ, has the largest instability region, followed in descending order by modes ðm ¼ 2; n ¼ 1Þ, ðm ¼ 1; n ¼ 1Þ, and ðm ¼ 3; n ¼ 1Þ, which has the smallest
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
0.50
Table 1 Comparison of unstable region for a panel having simply supported boundary condition under tensile and compressive loading ðb=R ¼ 0:5; L=R ¼ 2:0; N 0 ¼ 0:5N cr Þ Mode
0.45
Point of origin Tensile loading
0.03
0.04
0.05
Compressive loading
Ng et al. [15]
Present
Ng et al. [15]
Present
(1,1) (2,1) (3,1) (4,1)
0.7651 1.0722 1.5076 1.9976
0.7512 1.0661 1.4929 1.9973
0.6533 0.7182 0.9172 1.1958
0.6352 0.6956 0.8564 1.1371
(1,1) (2,1) (3,1) (4,1)
0.9977 1.3375 1.8313 2.4111
0.9756 1.3229 1.8140 2.4148
0.8854 0.9716 1.2024 1.5495
0.8839 0.9520 1.1409 1.4969
(1,1) (2,1) (3,1) (4,1)
1.2308 1.6063 2.1621 2.8345
1.2008 1.5833 2.1408 2.8405
1.1183 1.2347 1.5135 1.9419
1.0836 1.1921 1.4489 1.8929
0.40
mode 1 P
h=R
1161
0.35
mode 2 mode 3 mode 4
0.30
0.25
0.20 0.65
0.0
0.1
0.60
0.2
0.3 Ns / No
0.4
0.5
Fig. 3. Unstable regions for the first four modes of cylindrical panel with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 200Þ.
0.55
P
0.50
0.34
0.45
0.32 0.40 mode 1
mode 3
mode 2
mode 4
0.30
0.35
0.28 0.30
mode 1
0.26 P
mode 2
0.25 0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 2. Unstable regions for the first four modes of cylindrical panel with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
instability region. It is also noted that the region size associated with the fourth mode ðm ¼ 2; n ¼ 2Þ is approximately twice as wide as the principal region of the fundamental mode ðm ¼ 2; n ¼ 1Þ. Fig. 3 presents the first four instability regions of a simply supported shell panel 1 with thickness ratio h=R ¼ 200 . The first four modes are, respectively, modes ðm ¼ 1; n ¼ 2Þ, ðm ¼ 1; n ¼ 3Þ,
mode 3
0.24
mode 4 0.22 0.20 0.18 0.16 0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 4. Unstable regions for the first four modes of cylindrical panel with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 300Þ.
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1162
0.40
0.14 0.13
0.35
0.12 0.30
0.11 0.10 P
P
mode 1 0.25
mode 2
mode 1 mode 2 mode 3 mode 4
0.09
mode 3 mode 4 0.20
0.08 0.07
0.15
0.06 0.10
0.05 0.0
0.0
0.1
0.2 Ns / No
0.3
0.4
0.1
0.5
0.2 Ns / No
0.3
0.4
0.5
Fig. 7. Unstable regions for the first four modes of cylindrical panel with CFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 300Þ.
Fig. 5. Unstable regions for the first four modes of cylindrical panel with CFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
0.26 0.18
0.24 0.16
mode 1
mode 1 mode 2 mode 3 mode 4
0.22
mode 2 mode 3
0.14
0.20 P
P
mode 4 0.12
0.18 0.10
0.16 0.08
0.14 0.06 0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 6. Unstable regions for the first four modes of cylindrical panel with CFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 200Þ.
0.0
0.1
0.2 Ns / No
0.3
0.4
0.5
Fig. 8. Unstable regions for the first four modes of cylindrical panel with FFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1163
0.71
0.16 mode 1 mode 2 mode 3 mode 4
0.15
o
o
o
0 /90 /90 /0
o
o
o
o
o
0 /90 /0 /90
0.70
ϖ
0.14 0.69
P
0.13
0.12
0.68
0.11 0.67 0.0
0.1
0.2
0.3
0.4
0.5
Ns / No
0.10
0.09 0.0
0.1
0.2
0.3
0.4
0.5
Fig. 11. Unstable regions for (0 =90 =0 =90 and 0 =90 =90 =0 ) cylindrical panels with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
Ns / No 0.89
Fig. 9. Unstable regions for the first four modes of cylindrical panel with FFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 200Þ.
0.88 0.87
0.12
45o/-45o 45o/-45o/45o/-45o 45o/-45o/45o/-45o/45o/-45o
0.86 ϖ
0.13 mode 1 mode 2 mode 3 mode 4
0.85 0.84 0.83
0.11
P
0.82 0.0
0.1
0.2
0.3 Ns / No
0.4
0.5
0.10 Fig. 12. Unstable regions for simply supportive cylindrical panels with different layers ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ.
0.09
0.08
0.0
0.1
0.2
0.3
0.4
0.5
Ns / No Fig. 10. Unstable regions for the first four modes of cylindrical panel with FFFF boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 300Þ.
ðm ¼ 1; n ¼ 1Þ, and ðm ¼ 1; n ¼ 4Þ. It is clear that the first four modes in Fig. 3 are of different order as with those in Fig. 2 due to the change in panel thickness. It is also observed that the widths of the instability regions of panels are narrower as compared with the corresponding h=R ¼ 1 100 cases, and increase in the sequence of modes ðm ¼ 1; n ¼ 2Þ, ðm ¼ 1; n ¼ 3Þ, ðm ¼ 1; n ¼ 1Þ, andðm ¼ 1; n ¼ 4Þ. It is further noted that the sizes of the instability regions of the two lowest modes, ðm ¼ 1; n ¼ 1Þ and ðm ¼ 1; n ¼ 4Þ, have no significant difference. For a panel 1 with thickness ratio h=R ¼ 300 , corresponding results are
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
1164
frequency parameter o ¯ is defined as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rð1 n12 n21 Þ . o ¯ ¼ RP E 22
0.72 α = 15o α = 30o
(60)
o
α = 45
0.68
α = 60o o
ϖ
α = 75 0.64
0.60
0.56 0.0
0.1
0.2
0.3
0.4
0.5
Ns / No Fig. 13. Unstable regions for ða= a=a= aÞ cylindrical panels with SSSS boundary conditions ðy0 ¼ 60 ; L=R ¼ 2; R=h ¼ 100Þ. 1 shown in Fig. 4, and similar trends to that of the h=R ¼ 100 case are observed. However, the first four modes are now ðm ¼ 1; n ¼ 2Þ, ðm ¼ 1; n ¼ 3Þ, ðm ¼ 1; n ¼ 4Þ, and ðm ¼ 2; n ¼ 3Þ, respectively. Figs. 5–7 show the first four principal instability regions of the CFFF cylindrical panels with respective thickness 1 1 1 ratios of h=R ¼ 100 , h=R ¼ 200 , and h=R ¼ 300 . For the panel 1 having thickness ratio h=R ¼ 100 , the instability region of mode 4 is the widest and much larger than that of modes 1–3. The instability region for mode 2 is larger than that of mode 3 but less than that of mode 1. Further, it is noticed that the instability regions of modes 1 and 2 are overlapping as N s =N 0 exceeds 0.33. The first four unstable 1 regions of a panel with h=R ¼ 200 show different trends. From Fig. 6, it is observed that the width of the unstable regions reduces from mode 1 to 4, and some portions of the unstable regions of modes 1 and 2 are overlapping. This trend is also observed from Fig. 7 for the panel having 1 thickness ratio h=R ¼ 300 . The respective unstable regions of completely free panels 1 1 1 with thickness ratios h=R ¼ 100 , h=R ¼ 200 , and h=R ¼ 300 are depicted in Figs. 8–10. It is obvious that all three figures show similar general trends, and the widths of the instability regions for all four modes in each case have no significant differences. It is also found that certain portions of the unstable regions of all four modes mutually overlap 1 for the panels with thickness ratios h=R ¼ 200 and 1 h=R ¼ 300.
3.2. Laminated composite cylindrical panels The material properties of the laminated cylindrical panels in the present study are E 1 =E 2 ¼ 40, G 12 =E 2 ¼ 0:5, h=R ¼ 1=100, and n ¼ 0:25. The static axial compressive force is chosen to be N 0 ¼ 0:3N cr . The non-dimensional
For the simply supported panels with the symmetric and anti-symmetric cross-ply lamination schemes, 0 =90 =0 = 90 and 0 =90 =90 =0 , the unstable regions are plotted in Fig. 11. The modes for the two lamination schemes correspond to ðm ¼ 1; n ¼ 1Þ. It is observed that the two regions have comparatively same width, and the instability region for panel with 0 =90 =0 =90 scheme is at higher frequencies than that of panel with 0 =90 =90 =0 scheme. The effects of number of layers on the instability regions of laminated panels are demonstrated in Fig. 12. The twolayer, four-layer and six-layer panels, which have lamination schemes of 45 = 45 , 45 = 45 =45 = 45 , and 45 = 45 =45 = 45 =45 = 45 , respectively, are investigated. The corresponding mode in this case is ðm ¼ 1; n ¼ 1Þ. It is noted that all instability regions have higher frequencies and become narrower as the number of layers increases. It is also seen that the instability regions for four-layer and six-layer laminates begin to overlap as N s =N 0 exceeds 0.31. It is expected that the overlapping area of instability regions will become larger when the number of layers increases further. Fig. 13 illustrates the influence of lamination angles on the instability regions of laminates. The lamination scheme a= a=a= a is selected, five lamination angles, 15 , 30 , 45 , 60 , and 75 , are considered in this case. It is observed that, except mode ðm ¼ 1; n ¼ 1Þ corresponded by panel with lamination angle 60 , the mode ðm ¼ 1; n ¼ 2Þ, relates to the panels with other four lamination angles. It is also noted that the instability region of panel with angle 45 has the highest frequencies, followed by those of panels with lamination angles of 15 , 60 , 30 , and 75 . The instability region widens in the same lamination angle sequence.
4. Conclusions The dynamic stability of thin circular cylindrical panels under static and periodic axial forces has been investigated using the mesh-free kp-Ritz method. The mesh-free kernel particle estimate was successfully employed to approximate the 2D transverse displacement field. The Ritz minimization procedure was employed to obtain a system of Mathieu–Hill equations, and the principal instability regions were calculated via Bolotin’s first approximation. Numerical comparisons were made with those in literature, and the present methodology validated. The instability regions of composite laminated cylindrical panels having different lamination schemes and layers were investigated. It was found that the both the lamination scheme and number of layer influenced the positions and sizes of the instability regions.
ARTICLE IN PRESS K.M. Liew et al. / International Journal of Mechanical Sciences 49 (2007) 1156–1165
Acknowledgements The work described in this paper was fully supported by a grant from the Research Grants Councils of the Hong Kong Special Administrative Region, China (9041011 CityU 1160/05E). References [1] Bolotin VV. The dynamic stability of elastic systems. Holden-Day, San Francisco; 1964. [2] Vijayaraghavan A, Evan-Iwanowski RM. Parametric instability of circular cylindrical shells. TASME, Journal of Applied Mechanics 1967;31:985–90. [3] Yao JC. Nonlinear elastic buckling and parametric excitation of a cylinder under axial loads. TASME, Journal of Applied Mechanics 1965;29:109–15. [4] Birman V. Dynamic stability of unsymmetrically laminated rectangular plates. Mechanics Research Communications 1985;12:81–6. [5] Srinivasan RS, Chellapandi P. Dynamic stability of rectangular laminated composite plates. Computers and Structures 1986;24:233–8. [6] Moorthy J, Reddy JN. Parametric instability of laminated composite plates with transverse shear deformation. International Journal of Solids and Structures 1990;26:801–11. [7] Bert CW, Birman V. Dynamic stability of shear deformable antisymmetric angle-ply plates. International Journal of Solids and Structures 1987;23:1053–61. [8] Bert CW, Birman V. Parametric instability of thick, orthotropic, circular cylindrical shells. Acta Mechanica 1988;71:61–76. [9] Koval LR. Effects of longitudinal resonance on the parametric stability of an axially excited cylindrical shell. Journal of Acoustical Society of America 1974;55:91–7. [10] Argento A, Scott RA. Dynamic instability of laminated anisotropic circular cylindrical shells, part I: theoretical development. Journal of Sound and Vibration 1993;162:311–22. [11] Argento A, Scott RA. Dynamic instability of laminated anisotropic circular cylindrical shells, part II: numerical results. Journal of Sound and Vibration 1993;162:323–32. [12] Ng TY, Lam KY, Reddy JN. Dynamic stability of cross-ply laminated composite cylindrical shells, International. Journal of Mechanical Sciences 1998;40:805–23. [13] Lam KY, Ng TY. Dynamic stability of cylindrical shells subjected to conservative periodic axial loads using different shell theories. Journal of Sound and Vibration 1997;207:497–520.
1165
[14] Ganapathi M, Varadan TK, Balamurugan V. Dynamic instability of laminated composite curved panels using finite element method. Computers and Structures 1994;53:335–42. [15] Ng TY, Lam KY, Reddy JN. Dynamic stability of cylindrical panels with transverse shear effects. International Journal of Solids and Structures 1999;36:3483–96. [16] Kitipornchai S, Xiang Y, Wang CM, Liew KM. Buckling of thick skew plates. International Journal for Numerical Method in Engineering 1993;36:1299–310. [17] Liew KM, Xiang Y, Wang CM, Kitipornchai S. Flexural vibration of shear deformable circular and annular plates on ring supports. Computer Methods in Applied Mechanics and Engineering 1993;110:301–15. [18] Liew KM, Xiang Y, Kitipornchai S, Meek JL. Formulation of Mindlin-Engesser model for stiffened plate vibration. Computer Methods in Applied Mechanics and Engineering 1995;120: 339–53. [19] Cheung YK, Zhou D. Free vibrations of rectangular unsymmetrically laminated composite plates with internal line supports. Computers and Structures 2001;79:1923–32. [20] Cheung YK, Zhou D. Three-dimensional vibration analysis of cantilevered and completely free isosceles triangular plates. International Journal of Solids Structures 2002;39:673–87. [21] Zeng H, Bert CW. Free vibration analysis of discretely stiffened skew plates. International Journal of Stability and Dynamics 2001;1: 125–44. [22] Su GH, Xiang Y. A non-discrete approach for analysis of plates with multiple subdomains. Engineering Structures 2002;24: 563–75. [23] Lim CW, Liew KM. A pb-2 Ritz formulation for flexural vibration of shallow cylindrical shells of rectangular planform. Journal of Sound and Vibration 1994;173:343–75. [24] Liew KM, Feng ZC. Three-dimensional free vibration analysis of perforated super elliptical plates via the p-Ritz method. International Journal of Mechanical Sciences 2001;43:2613–30. [25] Chen JS, Pan C, Wu CT. Large deformation analysis of rubber based on a reproducing kernel particle method. Computational Mechanics 1997;19:211–27. [26] Reddy JN. An introduction to the finite element method. New York: McGraw-Hill; 1993. [27] Reddy JN. Applied functional analysis and variational methods in engineering. New York: McGraw-Hill; 1986; reprinted by Krieger, Melbourne, Florida, 1991. [28] Timoshenko SP, Gere JM. Theory of elastic stability. New York: McGraw-Hill; 1961.