Buckling analysis of simply supported composite laminates subjected to an in-Plane compression load by a novel mesh-ffree method

Vietnam Journal of Mechanics, VAST, Vol. 33, No. 2 (2011), pp. 65 – 78 BUCKLING ANALYSIS OF SIMPLY SUPPORTED COMPOSITE LAMINATES SUBJECTED TO AN IN-PLANE COMPRESSION LOAD BY A NOVEL MESH-FFREE METHOD Tinh Quoc Bui University of Siegen, Germany Abstract. Buckling analysis of composite laminates under an in-plane compression load based on the mesh - free Galerkin Kriging method is presented. The moving Kriging in- terpolation (MK) technique possessing the delta property is employed to construct the shape functions, and thus no special techniques for imposing the essential boundary con- ditions are required. The present formulation is based on the Krichhoff plate theory. The applicability, the accuracy and the effectiveness of the method are illustrated through a number of numerical examples. The results calculated by the proposed method are com- pared with those of existing reference solutions available in the literature and very good agreements are observed. It can be said that the proposed method can be considered as an alternative numerical technique in terms of meshfree methods. Keywords: Buckling, Composite laminates, mesh - free method. 1. INTRODUCTION Multi-layered angle-ply composite structures or laminates have been increasingly considered for a variety of engineering applications. The key advantages of the composite materials are due to their good characteristics with higher stiffness-to-weight, strength-to- weight ratios, etc. over the traditional materials. In modern advanced industries such as marine structures, automobile, naval, civil, construction sectors and aerospace, a thorough understanding of buckling phenomenon is necessary because of the reliability in design. Exact solutions for arbitrary geometries are very difficult to be obtained. Experiments are not only a time-consuming task but also an expensive procedure. Therefore, an approxi- mate solution based on the numerical computational approaches is unavoidable. Many methods have been introduced to buckling analysis of composite laminates. Typically, Reddy et al [1] proposed a plate bending element accounted for buckling and vibration of laminates. Huang and Li [2] developed the moving least square differential quadrature method (MLSDQ). Recently, Belinha et al [3] and Dai et al [4] used the element free Galerkin method while the mesh - free radial basis function method is developed by Liu et al [5]. The conforming radial point interpolation is studied by Liu et al [6]. For more information, readers can refer to various interesting reviews see e.g., [7-11]. 66 Tinh Quoc Bui The finite element method (FEM) is one of the most popular used tools in practice and has been developed for analyzing laminated composite plate problems [11]. Recently, a novel class of numerical methods named meshfree methods has been introduced [12- 16], where the problem domain is discretized by a set of nodes scattered regardless of the concepts of elements or meshing. Various meshfree approaches have been proposed and applied to analyses of composite laminates [17-22] involving the reproducing kernel particle method (RKPM), the meshless local Petrov-Galerkin (MLPG), the element-free Galerkin (EFG), the radial point interpolation method (RPIM) and so on. However, the lack of the delta property is awkward and encountered in most recent meshless methods and special techniques such as penalty method [12]; Lagrange multipliers [13]; coupling with the FEM [23], etc. are devoted to overcome that burdensome task. The moving Kriging interpolation method associated with the EFG method is first proposed by Gu [24] solving a simple problem of steady-state heat conduction. Further developments of the method can be found, respectively, for two-dimensional plane prob- lems [25, 26], shell structures [27], static deflections of thin plate [28], dynamic analysis of structures [29] and piezoelectric structures [30]. The objective of the paper is to extend the method to buckling analysis of laminated composite plates. A superior advantage of the method over the conventional approaches e.g., the moving least-square approximation (MLS) is that none of any special techniques needed for the enforcement of the essential boundary conditions are required. The structure of the paper is organized as follows. The meshfree formulations for buckling problem are briefly presented in the next section, in which the derivation of the MK shape functions, the governing equation and its discretiza- tion are involved. Section 3 is intended to present numerical results. We shall end with a conclusion in the last section. 2. MESH - FREE BUCKLING FORMULATION 2.1. Construction of shape functions The moving Kriging shape function is briefly presented in this section, see e.g., [24- 30] for more details. In the MK equation, a function u(x) in the sub-domain Ωx with Ωx ⊆ Ω is often approximated by u h(x), ∀x ∈ Ωx as uh(x) = [pT(x)A+ rT(x)B]u(x) (1) or uh(x) = n∑ I φI(x)uI (2) with uI = [u1 u2 ... un] T and n is the total number of the nodes in Ωx whereas φI(x) are the MK shape functions defined by φI(x) = m∑ j pj(x)AjI + n∑ k rk(x)BkI (3) The matrices A and B given in Eq. (3) are determined through A = (PTR−1P)−1PTR−1 (4) Buckling analysis of simply supported composite laminates subjected... 67 B = R−1(I−PA) (5) with I is an unit matrix, the vector p(x) is the polynomial with m basis functions p(x) = { p1(x) p2(x) · · · pm(x) }T (6) The matrix P is collected values of the polynomial basis functions and r(x) is given as r(x) = { R(x1, x) R(x2, x) · · · R(xn, x) } T (7) where R(xi, xj) is the correlation function between any pair of the n nodes xi and xj, it is denoted belong to the covariance of the field value u(x): R(xi, xj) = cov[u(xi)u(xj)] and R(xi, x) = cov[u(xi)u(x)]. Often, a Gaussian function with a correlation parameter θ is used to best fit the model R(xi, xj) = e −θr2ij (8) where rij = ‖xi − xj‖ and θ > 0 is defined as a correlation parameter which has a signif- icant influence on the solution. Thus, this parameter is considered in the numerical part. The quadratic basis function pT(x) = [ 1 x y x2 y2 xy ] is taken for all numerical computations. Because the second-order derivatives of the shape functions are required for analysis of thin plates, they are provided as follows φI,i(x) = m∑ j pj,i(x)AjI + n∑ k rk,i(x)BkI (9) φI,ii(x) = m∑ j pj,ii(x)AjI + n∑ k rk,ii(x)BkI (10) In mesh - free methods, the radius of influence domain used to determine the number of scattered nodes within an interpolated domain of interest is needed. Often, the following relation is employed to compute the size of support domain dm = αdc (11) where dc is a characteristic length regarding the nodal spacing close to the point of interest while α stands for a scaling factor. In practice, the size of the domain of influence must be large enough to sufficiently cover of the neighboring nodes. This implies that the scaling factor must be chosen some- how to ensure all necessary scattered nodes lying inside the domain of influences so that the problem can be converged. The optimal values of this scaling factor are depended on the problems of interest an optimal range from 2.0 to 4.0 is often taken [12, 13, 16]. On the other hand, the Gaussian correlation function in Eq. (8) is strongly sensitive to the correlation parameter whose value is found to be unrelated to any physical aspect of the problem. Deriving optimal values of the correlation parameter for all problems is very difficult. It varies from one to another problem and in theory no exact rules to get such a single optimal value for all problems. Hence, it is of interest to alternatively evaluate of the correlation parameter so that there should be existed an acceptable range on its magnitude to ensure consistency in the quality of the results. 68 Tinh Quoc Bui 2.2. Discrete governing equations Consider a laminate under a Cartesian coordinate system with the thickness t in the z direction consisting of nL ply layers as schematically represented in Fig. 1 showing the fiber direction of a layer denoted by ϕ. The governing equation of the general plate under buckling loads takes place as [10] ∂4w ∂x4 + 2 ∂4w ∂x2∂y2 + ∂4w ∂y4 = 1 D ( Nx ∂2w ∂x2 + 2Nxy ∂2w ∂x∂y +Ny ∂2w ∂y2 ) (12) where D denotes the matrix of constants relative to the laminated material property and the thickness. The laminated plate is subjected to in-plane compression forces (Nx, Ny, Nxy) that cause the instability are expressed as N = { Nx Ny Nxy } T with Nx = −N0;Ny = −µ1N0;Nxy = −µ2N0 (13) where N0 is a constant and µ1 and µ2 are possibly functions of coordinates [12]. Based on the classical plate theory, the deflection w(x) with x = {x, y}T is approximated using parameter of nodal deflection wI wh(x) = n∑ I φI(x)wI (14) where the φI(x) are the meshfree MK shape functions given in Eq. (3) above. The dis- placement fields are defined by u = { u v w }T = { −z ∂w ∂x −z ∂w ∂y w }T =Hw (15) with H = { −z ∂ ∂x − z ∂ ∂y 1 }T . Fig. 1. Geometry and its coordinate system of plies of a laminated composite plate showing a fiber orientation of ϕ at the top layer Buckling analysis of simply supported composite laminates subjected... 69 The strains and stresses of the plate are obtained by εp = { − ∂2w ∂x2 − ∂2w ∂y2 −2 ∂2w ∂x∂y }T = Lw (16) σp = { Mx My Mxy }T . (17) In Eq. (16), L = { − ∂2 ∂x2 − ∂2 ∂x2 − ∂2 ∂x∂y }T and in Eq. (17), Mx,My and Mxy are bending and twisting moments, respectively. The constitutive equation of the relationship between the strains and stresses can be thus expressed as σp = Dεp (18) where σp and εp are defined as pseudo-strains and pseudo-stresses. For a laminated com- posite plate based on the assumption of the classical theory gives [10, 12] D =   D11 D12 D16D12 D22 D26 D16 D26 D66   (19) DIJ = 1 3 nL∑ kL (Q¯IJ)(z 3 kL − z3kL−1); I, J = 1, 2, 6. (20) In Eq. (20), nL is defined by the number of layers of the laminated composite plate and the quantities Q¯IJ are determined as follows Q¯11 = Q11 cos 4 ϕ+ 2(Q12 + 2Q66) sin 2 ϕ cos2 ϕ+Q22 sin 4 ϕ (21) Q¯12 = (Q11 +Q22 − 4Q66) sin 2 ϕ cos2 ϕ+Q12(sin 4 ϕ+ cos4 ϕ) (22) Q¯16 = (Q11 −Q12 − 2Q66) sinϕ cos 3 ϕ+ (Q12 −Q22 + 2Q66) sin 3 ϕ cosϕ (23) Q¯22 = Q11 sin 4 ϕ+ 2(Q12 + 2Q66) sin 2 ϕ cos2 ϕ+Q22 cos 4 ϕ (24) Q¯26 = (Q11 −Q12 − 2Q66) sin 3 ϕ cosϕ+ (Q12 −Q22 + 2Q66) sinϕ cos 3 ϕ (25) Q¯66 = (Q11 +Q22 − 2Q12 − 2Q66) sin 2 ϕ cos2 ϕ+Q66(sin 4 ϕ+ cos4 ϕ) (26) Q11 = E1 1− ν12ν21 ;Q12 = ν12E2 1− ν12ν21 ;Q22 = E2 1− ν12ν21 ;Q66 = G12 (27) ν21E1 = ν12E2 (28) where E1, E2 are Young’s moduli parallel to and perpendicular to the fibers orientation, respectively, and G12 is the shear modulus while ν12, ν21 standing for the Poisson’s ratios. To derive the discretized governing equations for buckling analysis of laminated plate, the strain energy Πε describing the relation between the pseudo-stresses σp and pseudo-strains εp is given Πε = 1 2 ∫ Ω εTp σpdΩ (29) 70 Tinh Quoc Bui and the strain energy caused by in-plane forces of the plate is given by Πi = 1 2 ∫ Ω [ Nx ∂2w ∂x2 +Ny ∂2w ∂y2 + 2Nxy ∂2w ∂x∂y ] dΩ = 1 2 ∫ Ω εTpNdΩ (30) with the total energy functional Πν Πν = Πε +Πi = 1 2 ∫ Ω εTp σpdΩ+ 1 2 ∫ Ω εTpNdΩ (31) Using the standard Galerkin weak form, the discrete equation for buckling problems is obtained by substituting Eqs. (16) - (18) into Eq. (31), yielding∫ Ω ∂ ∂w (Lw)TD(Lw)dΩ+ ∫ Ω ∂ ∂w (Lw)TNdΩ = 0 (32) Substituting the deflection field w given in Eq. (14) into the variational form shown in Eq. (32) involving the effect of both the strain energies caused by bending and by in- plane forces, we finally obtain the discrete equation for buckling analysis of the laminated composite plate (K−N0K g)w = 0 (33) with K is the global stiffness matrix KIJ = ∫ Ω BTI DBJdΩ (34) BI = { − ∂2φI ∂x2 − ∂2φI ∂y2 −2 ∂2φI ∂x∂y }T (35) with N0 is the eigenvalue for an unitary compressive load or critical buckling loads which is needs to be determined, and Kg is the geometric stiffness matrix : K g IJ = ∫ Ω ∂φI ∂x ∂φJ ∂x + µ1 ∂φI ∂y ∂φJ ∂y + µ2( ∂φI ∂x ∂φJ ∂y + ∂φJ ∂x ∂φI ∂y )dΩ (36) 3. NUMERICAL EXAMPLES It must be noted here that only the fully simply supported boundary is employed and the fully clamped one is not taken into account because two rotations defined as unknown variable are not evolved in the approximation function. Therefore, it enables us to treat this problem with the present statement of the interpolation formulation. However, a Hermite-type technique developed in the RPIM [31, 32] might be applicable where both deflection and its derivatives are defined as variables field in the interpolation. The MK approach may have a similar manner which is thus needed such a development. However, the task would be more challenging and beyond the scope of the present work. Buckling analysis of simply supported composite laminates subjected... 71 3.1. Rectangular laminates A laminated rectangular plate made of E-glass/epoxy materials with its lengths a/b is considered. If not specified otherwise, the following material parameters are employed E1/E2 = 2.45, G12 = 0.48E2, the Poisson’s ratio ν12 = 0.23, the thickness t = 0.06m and the non-dimensional buckling load k = N0b 2/pi2D0 with D0 = E1t 3/(12(1− ν12ν21)). The in-plane compressive load Nx = −N0 applied in the x direction is considered only. 3.1.1. Effect of the scaling and correlation parameters For convenience in evaluating the results with other numerical approaches, the plate is set to be a square one a = b = 10m. As mentioned previously, the scaling factor α and the correlation parameter θ have a significant effect on the solution. They are therefore considered first in order to find out an acceptable range to be used. A regular set of 13×13 scattered nodes is taken. For the scaling factor analysis, a correlation parameter of 3.5 is typically specified and unchanged while the α is varied within a wide range 2.0 ≤ α ≤ 6.0. The results of the dimensionless critical buckling loads corresponding to each specified value of the scaling factor are presented in Table 1. They are compared with the results obtained by the EFG. Based on the achieved results, a value of 2.5 ≤ α ≤ 3.5must be used in practice. For the correlation parameter analysis, Table 2 presents the influence of the parameter on the buckling loads. The correlation parameter is varied in a specified wide range, i.e., 1 ≤ θ ≤ 10 whereas a scaling factor of 3 is fixed. The results of the dimensionless critical buckling loads corresponding to each specified value of the correlation parameter are compared with that of the EFG [12]. It is found that a value of 3 ≤ θ ≤ 4 must be chosen in practice so that a reasonable critical buckling load can be obtained. Table 1. Evaluation of the scaling factor on the critical buckling loads with various angle-ply orientations for the laminated composite square plate Angle-ply The scaling factor EFG [12] 2.0 2.5 2.8 3 3.5 4.0 6.0 (00, 00, 00) 3.648 2.361 2.396 2.391 2.553 2.960 6.645 2.39 (150,−150, 150) 3.623 2.431 2.442 2.448 2.580 2.945 6.470 2.45 (300,−300, 300) 3.627 2.530 2.540 2.578 2.655 2.977 6.199 2.57 (450,−450, 450) 3.642 2.595 2.587 2.649 2.671 2.964 5.977 2.64 (00, 900, 00) 3.671 2.348 2.380 2.394 2.547 2.939 6.632 2.39 3.1.2. Effect of the length-to-width ratio The angle-ply of three layers arranged as (00, 900, 00) is again used associated with various regular patterns of discretized nodes. The length-to-width ratio a/b is varied by specifying various typical magnitudes and its corresponding calculated results for the crit- ical buckling loads are presented in Table 3 and Fig. 2. It is found that the dimensionless buckling loads are decreased once the length-to-width ratio a/b decreasing from 2.5 to 0.8. Furthermore, Fig. 3 provides the first buckling mode shape for different length-to-width ratios a/b = 1.0; 1.5; 2.0 and 2.5, respectively. 72 Tinh Quoc Bui Table 2. Evaluation of the correlation parameter on the critical buckling loads with various angle-ply orientations for the laminated composite square plate Angle-ply The correlation parameter EFG [12] 1 2 3 3.5 4 5 10 (00, 00, 00) 2.172 2.371 2.367 2.391 2.437 2.678 3.522 2.39 (150,−150, 150) 2.212 2.389 2.411 2.452 2.562 2.811 3.568 2.45 (300,−300, 300) 2.301 2.455 2.495 2.578 2.680 2.930 3.777 2.57 (450,−450, 450) 2.340 2.476 2.588 2.667 2.735 2.966 3.864 2.64 (00, 900, 00) 2.277 2.307 2.389 2.391 2.455 2.549 3.215 2.39 Table 3. Dimensionless buckling loads versus the length-to-width ratio a/b for the laminated composite plates a/b Nodes 5× 5 7× 7 9× 9 11× 11 13× 13 15× 15 2.5 15.0816 11.1746 7.1746 4.9746 4.8746 4.9859 2.0 16.3567 11.0228 6.9570 3.9889 3.8364 3.8409 1.5 23.9234 10.3261 6.5171 3.5991 3.1549 3.2608 1.2 27.9351 9.6677 4.4861 2.9520 2.7097 2.7139 1.0 26.9469 8.1661 4.0127 2.5299 2.5857 2.5894 0.8 16.4483 5.4583 3.6498 2.5248 2.5868 2.5274 Fig. 2. Convergence of dimensionless buckling loads versus the length-to-width ratio a/b through the number of scattered nodes for the laminated composite plates Buckling analysis of simply supported composite laminates subjected... 73 Fig. 3. The first buckling mode of rectangular laminated composite plate at var- ious length-to-width a/b = 1 (top-left), a/b = 1.5(top-right), a/b = 2.0 (bottom- left) and a/b = 2.5 (bottom-right) 3.1.3. Effect of the number of layers To further validate the code, a study of the effect of the number of layers of compos- ite laminates on the buckling loads is also examined. The dimensionless critical buckling loads are calculated for both regular and irregular nodes for four and five-ply laminates. The results are presented in Table 4, we found that the variation of the number of layers has a tiny effect on the critical buckling loads. Further, Fig. 4 also provides the first twenty buckling mode shapes of the laminated square plate with angle-ply of (00, 00, 00) to get a better observation. Table 4. Effect of number of layers on the dimensionless critical buckling loads for the laminated square plate with various angle-ply orientations Angle-ply Regular Irregular Four layers (00, 00, 00, 00) 2.4113 2.4139 (300,−300, 300,−300) 2.6240 2.6301 (450,−450, 450,−450) 2.6919 2.7032 (00, 900, 00, 900) 2.4103 2.4131 Five layers Regular Irregular (00, 00, 00, 00, 00) 2.4113 2.4139 (300,−300, 300,−300, 300) 2.5981 2.6004 (450,−450, 450,−450, 450) 2.6578 2.6691 (00, 900, 00, 900, 00) 2.4113 2.4139 74 Tinh Quoc Bui Fig. 4. The first twenty buckling mode shapes of the laminated composite square plate with angle-ply of (00, 00, 00) 3.1.4. Effect of the modulus ratio E1/E2 In this subsection, a study of the influence of modulus ratioE1/E2 on the critical buckling loads is considered. The material properties of the laminated plate are as follows: a = b = 10m in length, the thickness a/t = 10, other ratios concerning the elastic constants as E1/E2 varied, G12 = 0.6E2 and the Poisson ratio ν12 = 0.25. The dimensionless critical buckling factor is calculated by k = N0a 2/t3E2. The modulus ratio E1/E2 is varied by several specified values, respectively, to that the dimensionless buckling loads coefficients are calculated correspondingly. The computed results for both regular and irregular nodes are listed in Table 5 and Fig. 5, respectively. These results are compared with those studied by Liu et al. [5] and Phan and Reddy [33] using the same classical plate theory (CLPT). An excellent agreement is obtained. The buckling loads are increased once the modulus ratios E1/E2 increasing. 3.2. Laminated composite square plate with a hole of complicated shape Another square plate with a hole of complicated shape is also tested to illustrate the applicability of the method to arbitrary geometries. The geometry is depicted in Fig. 6 including its nodal distribution. Two irregular patterns of 134 and 506 scattered nodes are taken. The material parameters are the same as the previous examples. For reference Buckling analysis of simply supported composite laminates subjected... 75 Table 5. Comparison of the dimensionless buckling loads among approaches af- fected by the modulus ratio for the laminated squared plate E1/E2 CLPT Liu et al [5] CLPT Phan et al [33] Present Regular Irregular 3 5.761 5.7358 5.7740 5.7792 10 11.576 11.492 11.5382 11.5753 20 20.127 19.712 20.1985 20.2747 30 28.232 27.936 28.8606 28.8878 40 36.367 36.160 36.5231 36.5798 Fig. 5. Comparison of the dimensionless critical buckling loads varied by the mod- ulus ratios Fig. 6. A laminated composite square plate with a hole of complicated shape with 506 nodes and its geometry 76 Tinh Quoc Bui solutions, an extra task calculated by the MLS-based EFG using the Lagrange multipliers is also coded in the same manner to that the present method does as well as a fine FEM solution solved by ANSYS software with 15288 DOFs. The results of the dimensionless buckling loads are presented in Table 6 and Fig. 7 with various orientations of three-ply in a comparison with the EFG’s reults. A very good agreement to each other is obtained with α = 3 and θ = 2, 3 and 3.5 are used, respectively. Table 6. Comparison of the dimensionless buckling loads for the laminated square plate with a hole of complicated shape between the present and the EFG methods Angle-ply 134 nodes 506 nodes FEM (ANSYS) EFG Present EFG Present θ = 2.0 3.0 3.5 (00, 00, 00) 1.2048 1.2132 1.2048 1.2322 1.2210 1.2232 1.2107345 (150,−150, 150) 1.3358 1.3169 1.2986 1.3112 1.3162 1.3172 1.3091655 (300,−300, 300) 1.5938 1.5861 1.5524 1.5719 1.5754 1.5763 1.5763909 (450,−450, 450) 1.7176 1.6838 1.6733 1.6903 1.6856 1.6878 1.6742200 (00, 900, 00) 1.2048 1.2082 1.2048 1.2801 1.2427 1.2433 1.2121556 Fig. 7. Comparison of the dimensionless buckling loads for the laminated square plate with a hole of complicated shape among the present, the EFG methods and the FEM 4. CONCLUSIONS A detailed buckling analysis of laminated composite plates under an in-plane com- pression load using the mesh - free Galerkin Kriging method is presented. The applicability and the accuracy of the method are demonstrated through a number of solved numerical examples comparing the results with existing solutions. Very good agreements have been observed. The effect of various ratios aspects such as the fiber orientations, the modu- lus, the scaling and correlation parameters, length-to-width and number of layers on the Buckling analysis of simply supported composite laminates subjected... 77 buckling loads is examined. It is seen that big advantage when using this method is that the difficulty of the enforcement of the boundary condition can be avoided completely, whereas the optimal choice of the correlation parameter for all problems is known as its disadvantage. As a consequence, the method is marginal improvement and adequately accurate and its applications to other complex problems are of course promising. REFERENCES [1] Reddy J. N., Khdeir A. A., Buckling and vibration of laminated plates using various plate theories, AIAA Journal, 27, (1989), 1808-1817. [2] Huang Y. Q., Li Q. S., Bending and buckling analysis of antisymmetric laminates using the moving least square differential quadrature method, Computer Methods in Applied Mechanics and Engineering, 193, (2004), 3471-3492. [3] Belinha J., Dinis LMJS, Analysis of plates and laminates using the element free Galerkin method, Computer and Structures, 84, (2006), 1547-1559. [4] Dai K. Y., Liu G. R., Lim K.M., Chen X.L., A mesh-free method for static and free vibration analysis of shear deformable laminated composite plates, Journal of Sound and Vibration, 269, (2004), 633-652. [5] Liu L., Chua L. P., Ghista D.N., Mesh-free radial basis function method for static, free vi- bration and buckling analysis of shear deformable composite laminates, Composits Structures, 78, (2007), 58-69. [6] Liu G. R., Zhao X., Dai K. Y., Zhong Z. H., Li G. Y., Han X., Static and free vibration analysis of laminated composite plates using the conforming radial point interpolation method, Composite Science and Technology, 68, (2008), 354-366. [7] Kapania R. K., Raciti S., Recent advances in analysis of laminated beams and plates, part I: shear effects and buckling, AIAA Journal, 27, (1989), 923-934. [8] Zhang Y. X., Yang C. H., Recent developments in finite element analysis for laminated com- posite plates, Composite Structures, 88, (2009), 147-157. [9] Reddy J. N., Miravete A., Practical analysis of composite laminates, Boca Raton (FL): CRC Press, (1995). [10] Reddy J. N., Mechanics of laminated composite plates and shells: theory and analysis, CRC Press, (1996). [11] Ochoa O. O., Reddy J. N., Finite element analysis of composite laminates, Kluwer Academic Publishers, USA, (1992). [12] Liu G. R., Meshfree methods: moving beyond the finite element method, CRC Press: U.S.A, (2003). [13] Belytschko T., Lu Y. Y., Gu L., Element free Galerkin method, International Journal for Numerical Methods and Engineering, 37, (1994), 229-256. [14] Atluri S. N., Zhu T., A new meshless local Petrov-Galerkin (MLPG) approach in computa- tional mechanics, Computational Mechanics, 22, (1998), 117-127. [15] Liu W. K., Jun S., Zhang Y. F., Reproducing kernel particle method, International Journal for Numerical Methods in Fluids, 20, (1995), 1081-1106. [16] Li S., Liw W. K., Meshfree particle method, Springer Berlin-Heidelberg New York, (2004). [17] Wang J., Liew K. M., Tan M. J., Rajendran S., Analysis of rectangular laminated composite plates via FSDT meshless method, International Journal of Mechanical Sciences, 44, (2002), 1275-1293. 78 Tinh Quoc Bui [18] Xiao J. R., Gilhooley D. F., Batra R. C., Gillespie J. W., Mccarthy M. A., Analysis of thick composite laminates using a higher-order shear and normal deformable plate theory (HOS- NDPT) and a meshless method, Composite Part B-Engineering, 39, (2008), 414-427. [19] Liew K. M., Lim H. K., TanM. K., He X. Q., Analysis of laminated composite beams and plates with piezoelectric patches using the element-free Galerkin method, Computational Mechanics, 29, (2002), 486-497. [20] Amirani M. C., Khalili S. M. R., Nemati N., Free vibration analysis of sandwich beam with FG core using the element free Galerkin method, Composite Structures, 90, (2009), 373-379. [21] Ferreira A. J. M., Roque C. M. C., Jorge R. M. N., Free vibration analysis of symmetric laminated composite plates by FSDT and radial basis functions, Computer Methods in Applied Mechanics and Engineering, 194, (2005), 4265-4278. [22] Xiang S., Wang K. M., Ai Y. T., Sha Y. D., Shi H., Analysis of isotropic, sandwich and laminated plates by a meshless method and various shear deformation theories, Composite Structures, 91, (2009), 31-37. [23] Belyschko T., Organ D., Krongauz Y., A coupled finite element - element frer Galerkinmethod, Computational Mechanics, 17, (1995), 186-195. [24] Gu L., Moving Kriging interpolation and element free Galerkin method, International Journal for Numerical Methods and Engineering, 56, (2003), 1-11. [25] Tongsuk P., Kanok-Nukulchai W., On the parametric refinement of moving Kriging interpola- tions for element free Galerkin method, Computational Mechanics WCCM VI in conjunction with APCOM’04, 2004; Sept. 5-10, Beijing, China. [26] Tongsuk P., Kanok-Nukulchai W., Further investigation of element free Galerkin method using moving Kriging interpolation, International Journal of Computational Methods, 1, (2004), 1-21. [27] Sayakoummane V., Kanok-Nukulchai W., A meshless analysis of shells based on moving Krig- ing interpolation, International Journal of Computational Methods, 4, (2007), 543-565. [28] Bui Q. T., Nguyen N. T., Nguyen-Dang H., A moving Kriging interpolation- based meshless method for numerical simulation of Kirchhoff plate problems, International Journal for Numerical Methods and Engineering, 77, (2009), 1371-1395. [29] Bui Q. T., Nguyen N. M., Zhang Ch., A moving Kriging interpolation-based element-free Galerkin method for structural dynamic analysis, Computer Methods in Applied Mechanics and Engineering, 200 (2011) 1345 - 1366. [30] Bui Q. T., Nguyen N. M., Zhang Ch., Pham D. A. K., An efficient meshfree method for analysis of two-dimensional piezoelectric structures, Smart Materials and Structures, 20(6) (2011). [31] Liu Y., Hon Y. X., Liew K. M., A meshfree Hermite-type radial point interpolation method for Kirchhoff plate problems, International Journal for Numerical Methods and Engineering, 66, (2006), 1153-1178. [32] Cui X., Liu R. G., Li G., A smoothed Hermit radial point interpolation method for thin plate analysis, Archive Applied Mechanics, (2009), DOI 10.1007/s00419-009-0392-0 [33] Phan N. D., Reddy J. N., Analysis of laminated composite plates using a higher-order shear deformation theory, International Journal for Numerical Methods and Engineering, 21, (1985), 2201-2219. Received September 14, 2010