- Email: [email protected]
Bending and vibration responses of laminated composite plates using an edge-based smoothing technique
Engineering Analysis with Boundary Elements 35 (2011) 818–826
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Bending and vibration responses of laminated composite plates using an edge-based smoothing technique X.Y. Cui a, G.R. Liu b, G.Y. Li a,n a b
State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, PR China School of Aerospace Systems, University of Cincinnati, Cincinnati, OH 45221-0070, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 November 2009 Accepted 16 December 2010 Available online 25 February 2011
In this paper, bending and vibration analysis of laminated composite plates is carried out using a novel triangular composite plate element based on an edge-based smoothing technique. The present formulation is based on the first-order shear deformation theory, and the discrete shear gap (DSG) method is employed to mitigate the shear locking. The smoothed Galerkin weak form is adopted to obtain the discretized system equations, and edge-based smoothing domains are used for the numerical integration to improve the accuracy and the convergence rate of the method. The present formulation is coded and used to solve various example problems of bending and free vibration of laminated composite plates. It is found that the present method can provide excellent results with a wide range of thickness and is free of shear locking. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Numerical methods Vibration FEM Techniques-FEM Laminated composite plates
1. Introduction The bending and vibration analysis of laminated composite plates has become more and more important in the past several decades, due to rapidly increased application of composite materials in structure systems in civil, mechanical and aerospace engineering. The search for simple, efficient and inexpensive numerical methods applicable to both thin and thick laminated composite plates has been receiving much interest. Various laminated plate theories have been introduced for analysis of thin to thick laminated plates such as the classical laminated plate theory (CLPT), the first-order shear deformation theory (FSDT) and the higher-order shear deformation theory (HSDT) [1]. Using these theories, some existing numerical methods have been adopted to compute bending deformation and natural frequencies of composite laminates. Reddy [2] developed the free vibration analysis of antisymmetric, angle-ply laminated plates including transverse shear deformation using the finite element method. A 3-node higher-order triangular plate element TRIPLT was formulated by Lakshminarayana and Murthy [3]. This element has 15 degrees of freedom (DOF) containing displacements and rotations along with their first derivatives per node. Noor and Burton [4] gave an assessment of computational models for multilayered anisotropic plates. Dai et al. [5] studied the static and free vibration analysis based on HSDT using element free Galerkin
n
Corresponding author. Tel.: +86 731 8821717; fax: + 86 731 8822051. E-mail address: [email protected] (G.Y. Li).
0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.01.007
(EFG) method. Setoodeh and Karami [6] investigated the characteristics of anisotropic thick laminated composite plates using a 3-D layer-wise FEM. Lee and Han [7] developed a 9-node assumed strain element for free and forced vibration analysis of laminated composite plates and shells. Although large numbers of laminated plate elements have been developed, most of them have a complex form. Displacement-based triangular elements have a very simple form as they only consider the displacement and rotation degrees of freedom at nodes and without additional parameters. However, these elements generally have a low accuracy because of the inherent characteristic of overstiff phenomenon from the standard weak formulation. In order to ‘‘soften’’ the system and allow the use of discontinuous shape functions, Liu [8] proposed a general formulation of smoothed Galerkin weak form using the strain smoothing technique [9]. Efficient numerical methods have been developed based on both FEM and meshfree settings using different smoothing domains that can be node-based cell-based can edge-based. Using the point interpolation method and the node-based smoothing operations, a node-based smoothed point interpolation method (NS-PIM) were developed using polynomial basis functions (NS-PIM [10]) and radial basis functions (LC-RPIM [11]). Liu and Zhang [12] proved that the NS-PIM is variationally consistent, and can provide much better stress results. More importantly it can provide upper bound solution in energy norm. Liu et al. [13] has found that the LC-RPIM needs a high order interpolation to obtain satisfactory solutions for laminated composite plates. However, the node-based smoothed models are ‘‘overly-soft’’, which leads to temporal instability [14–16]. The instability can be clearly shown as spurious non-zero energy modes
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
in free vibration analyses. Using the FEM shape functions and strain smoothing operations based on cells that are obtained by dividing the quadrilateral elements (Q4), a smoothed finite element method (SFEM) was first formulated by Liu and co-workers [17]. Liu et al. [18] gave detailed theoretical aspects including stability, bound property and convergence about SFEM, and Cui et al. [19] extended the SFEM for linear and geometrically nonlinear analysis of plates and shells. The SFEM idea was also applied to the MITC4 element by Nguyen-Xuan et al. [20] for linear and static analysis and free vibration analysis [21] of plates. The formulation of SFEM is found much flexible and robust than the FEM, and was further extended to general n-sided polygonal elements [22]. However, the SFEM does not give any advantage, when triangular elements (T3) are used, because the SFEM-T3 is in fact exactly the same as the FEM-T3. Recently, an edge-based smoothed finite element (ES-FEM) is proposed by constructing edge-based smoothing domain using the smoothed Galerkin weak form. The ES-FEM exhibit superconvergence properties in both energy norm and displacement norm and can produce ultra-accurate solutions in 2D solid mechanics problems [23]. In this work, the ES-FEM is further extended to solve bending and vibration problems of laminated composite plates using simple 3-node triangular elements that can be generated automatically. The present formulation is based on the first-order shear deformation theory, and shear locking effect is suppressed using the discrete shear gap method [24]. In the present method, the smoothing domains associated with the edges of the triangles are formed based on the triangular element mesh. The compatible strains in the smoothing domains are smoothed using the simple rectangle rule of integration, and the system stiffness matrix is obtained by numerical integration based on smoothing domains. The present method has a very simple form and without using any internal parameters in the process of building the stiffness matrix. Using the present ES-FEM method, excellent results have been obtained for both bending and vibration problems of laminated composite plates.
819
of the transverse normals, to the undeformed mid-plane, in the xz and yz planes, respectively. The strain vector e can be written in terms of the mid-plane deformations: ( ) 0 em zeb e ¼ fexx , eyy , gxy , gxz , gyz gT ¼ þ þ ð2Þ e 0 0 s where em is the membrane strain, eb is the bending strain and es is the shear strain, which can be given by 9 9 8 8 @y @u0 x 8 9 > > > > > > > > @x @x > > > > 0 = = < < < @w = @yy @v0 @x yx , eb ¼ , es ¼ @w0 ð3Þ em ¼ @y @y > > > > : @y yy ; > > > > > > ; ; : @u0 þ @v0 > : @yx þ @yy > @y
@x
@y
@x
The constitutive relationship of laminated composite plates can be expressed as
r^ ¼ D^ e^
ð4Þ
where 2
r^ ¼ f r^ Tm r^ Tb r^ Ts gT , e^ ¼ f eTm eTb e
T T s g ,
3
A
B
0
^ ¼6 D 4B 0
D
7 05
0
S
ð5Þ
^ m ¼ fNx ,Ny ,Nxy gT is the membrane force vector in the midwhere r ^ b ¼ fMx ,My ,Mxy gT is the bending moment vector and plane, r r^ s ¼ fQx , Qy gT is the transverse shear force vector. The membrane stiffness (A), the bending stiffness (D), the bending–extensional coupling stiffness (B) and the transverse shear stiffness (S) are defined as 2 3 2 3 A11 A12 A16 B11 B12 B16 6A 7 6 7 A ¼ 4 12 A22 A26 5, B ¼ 4 B12 B22 B26 5 A16 A26 A66 B16 B26 B66 2 3 " # D11 D12 D16 pffiffiffiffiffiffiffiffiffiffiffi k1 S44 k1 k2 S45 6D 7 D D p ffiffiffiffiffiffiffiffiffiffiffi ð6Þ S¼ D ¼ 4 12 22 26 5, k1 k2 S45 k2 S55 D16 D26 D66
2. Basic formulations in which k1 and k2 are shear correction factors. The stiffness in Eq. (6) can be defined as
2.1. Geometry and kinematics of laminated composite plate A typical laminated composite plate with thickness t and Nl layers is shown in Fig. 1. Based on the first-order shear deformation laminated plate theory [1], the displacements field can be expressed by uðx,y,zÞ ¼ u0 ðx,yÞ þzyx ðx,yÞ vðx,y,zÞ ¼ v0 ðx,yÞ þ zyy ðx,yÞ wðx,y,zÞ ¼ w0 ðx,yÞ
ð1Þ
in which u0, v0 and w0 are the displacements of the mid-plane of the plate in the x, y and z directions, yx and yy denote the rotations
Aij ¼ Dij ¼
Nl X
ðkÞ
Q ij ðtk tk þ 1 Þ,
k¼1 Nl X
1 ðkÞ Q ðt 3 t 3 Þ, 3 k ¼ 1 ij k k þ 1
Bij ¼
Nl 1X ðkÞ Q ðt 2 t 2 Þ 2 k ¼ 1 ij k k þ 1
Sij ¼
Nl X
ðkÞ
Q ij ðtk tk þ 1 Þ
ð7Þ
k¼1
where Aij, Bij and Dij are defined for i, j ¼1, 2, 6, whereas Sij is defined for i, j ¼4, 5. tk and tk + 1 denote the distances from the plate reference mid-plane to the outer and inner surfaces of the kth layer, respectively, as shown in Fig. 2. Nl is the total number of ðkÞ
layers in the laminated plate, and Q ij (i, j¼1, 2, 4, 5, 6) are related to the engineering constants for the kth layer and can refer to the work by Reddy [1].
t
2.2. Edge-based strain smoothing technique
Fig. 1. A typical laminated composite plate with thickness t and Nl layers and its coordinate system.
In the present method, the domain discretization is based on triangular elements which are same as those used in standard FEM. The problem domain O is divided into Ne triangular elements with a total of Nedge edges. The smoothing domain for each edge is formed by sequentially connecting two end points of the edge and centroids of its surrounding triangles to perform the smoothing operation, such that O ¼ O1 [ O2 [ [ ONedge and Oi\Oj ¼| (iaj, i¼1, y, Nedge, j¼ 1, y, Nedge). A strain smoothing
820
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
Fig. 3. The integration domain Ok of edge k is assembled by two sub-domains, Ok1 and Ok2. The sub-domain Ok1 is from the element e1, and the sub-domain Ok2 is form the element e2, xCk1 and xCk2 are the geometric centers of the sub-domains. Fig. 2. Cross-sectional view of the laminated plate.
technique is used over each smoothing domains and the integration is performed based on the smoothing domains associated with the edges, instead of based on the elements as in the standard FEM. The displacement at any point in the problem domain is interpolated using the element shape functions, which are created as same as in the linear standard FEM. The generalized mid-plane displacements vector u¼{u0, v0, w0, yx, yy}T consisting of both displacements and rotations can be then approximated as uðxÞ ¼
3 X
Ni ðxÞdi
ð8Þ
i¼1
where di ¼ {u0i, v0i, wi, yxi, yyi}T is the displacement of the node i, and Ni is the diagonal matrix of shape functions n
zfflfflfflfflfflfflffl}|fflfflfflfflfflfflffl{ Ni ¼ diagfg fNi ,. . .,Ni g
ð9Þ
in which n ¼5 is the degree of freedom in one node and Ni is the linear shape function. Introducing the strain smoothing operation [8], the strain in the smoothing domain Ok is given by Z 1 ek ¼ e^ ðxÞdO ð10Þ Ak Ok where e^ ðxÞ is the generalized compatible strain given in Eq. (4) and Ak is the area of the smoothing domain Ok. Since linear shape functions are used in the present method, the compatible strain e^ ðxÞ is a constant in each smoothing subdomain OkI. Therefore, the integration in Eq. (10) can be performed exactly using the simple rectangle rule of integration:
ek ¼
1 X AkI e^ ðxCkI Þ Ak
ð11Þ
where xCkI is the centroid of the smoothing sub-domain OkI, and AkI is the area of the smoothing sub-domain OkI. It is clear that all one needs to do is to sample the compatible strains at the geometric center of the smoothing sub-domain, which can be done very efficiently. As shown in Fig. 3, for interior edges, the integration domain Ok of edge k is formed by two sub-domains, Ok1 and Ok2. The subdomains Ok1 and Ok2 are from the element e1 and the element e2, respectively. The nodes of elements e1 and e2 are i, j, m and j, n, m, respectively. Substituting Eq. (8) into Eq. (3) and using the strain smoothing technique, the smoothed membrane strain can be rewritten as 1 emk ¼ ½Ak1 em ðxCk1 Þ þ Ak2 em ðxCk1 Þ Ak A A ðx Þd þ k2 Bk2 ðx Þd ¼ Bmk dk ¼ k1 Bk1 Ak m Ck1 k1 Ak m Ck2 k2
where dk1 and dk2 are the nodal displacement vectors of the element e1 and the element e2, respectively, dk are the displacements of the nodes associated with edge k, and BkI m ðxCkI Þ are the compatible strain matrices of the smoothing sub-domain OkI (I ¼1, 2), which are given as h i kI BkI BkI ðI ¼ 1,2Þ BkI m2 m3 , m ðxCk1 Þ ¼ Bm1 2 @N ðx Þ 3 i CkI 0 0 0 0 @x 6 7 @Ni ðxCkI Þ 6 0 0 0 07 BkI ð13Þ 7, ði ¼ 1,2,3Þ @y mi ¼ 6 4 5 @Ni ðxCkI Þ @Ni ðxCkI Þ 0 0 0 @y @x From Eq. (12), the smoothed membrane strain matrix Bmk can be written as Bmk ¼
Ak1 k1 Ak2 k2 B þ B Ak m Ak m
ð14Þ
Note that the sign ‘+ ’ denotes assembly but not summation here. Similarly, we can obtain the bending strain matrix Bbk as follows: Bbk ¼
Ak1 k1 Ak2 k2 B þ B Ak b Ak b
ð15Þ
in which h kI BkI b ðxCk1 Þ ¼ Bb1 2 0 0 0 6 60 0 0 kI Bbi ¼ 6 4 0 0 0
BkI b2
i BkI b3 ,
@Ni ðxCkI Þ @x
0 @Ni ðxCkI Þ @y
0
ðI ¼ 1,2Þ 3 7
@Ni ðxCkI Þ 7 7, @y @Ni ðxCkI Þ @x
ði ¼ 1,2,3Þ
5
The smoothed bending strain can be given as follows:
ebk ¼ Bbk dk
ð17Þ
In order to eliminate the shear locking, the ‘‘discrete shear gap’’ (DSG) method [24] is used. In each triangular element, the shear strain here can be computed by
gxz ¼
3 X @Ni i¼1
gyz ¼
@x
3 X @Ni i¼1
@y
3 X @Ni
Dwxi þ
i¼1
@x
3 X @Ni
Dwxi þ
i¼1
@y
Dwyi Dwyi
ð18Þ
where Dwxi and Dwyi are the discrete shear gaps at the node i, they are given by
Dwx1 ¼ Dwx3 ¼ Dwy1 ¼ Dwy2 ¼ 0 ð12Þ
ð16Þ
1 2
1 2
Dwx2 ¼ ðw2 w1 Þ þ aðyx1 þ yx2 Þ þ bðyy1 þ yy2 Þ
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
1 2
1 2
Dwy3 ¼ ðw3 w1 Þ þ cðyx1 þ yx3 Þ þ dðyy1 þ yy3 Þ
821
ð19Þ
in which a ¼ x2 x1 ,
b ¼ y2 y1
c ¼ x3 x1 ,
d ¼ y3 y1
ð20Þ
From Eqs. (18)–(20), the shear strain es in each element can be written as n oT es ¼ gxz gyz ¼ Bs de ð21Þ where de is the nodal displacement of the triangular element, Bs is the shear strain matrix given as Bs ¼ Bs1 Bs2 Bs3 " # 1 0 0 bd Ae 0 Bs1 ¼ 2Ae 0 0 ca 0 Ae " # ad bd 1 0 0 d 2 2 Bs2 ¼ 2Ae 0 0 c ac bc 2 2 " # bc bd 0 0 b 1 2 2 Bs3 ¼ ð22Þ ac ad 2Ae 0 0 a 2 2
Fig. 4. The integration domain Ok of edge k is assembled by only one domain.
matrix form: Kdf ¼ 0
ð28Þ
where f is the force vector defined as Z f ¼ NðxÞf~ dO
ð29Þ
O
and the smoothed stiffness matrix K are expressed as Nedge
For the smoothing domain Ok assembled by sub-domains Ok1 and Ok2, The shear strain matrices of the correlative triangular k2 element Bk1 s and Bs are computed from Eq. (22). The smoothed k shear strain es of the domain Ok can be written as
esk ¼ Bsk dk
ð23Þ
where Bsk is assembled by Bsk ¼
Ak1 k1 Ak2 k2 B þ B Ak s Ak s
Bk1 s
and
Bk2 s
as follows:
T
Bbk
T
Bsk
iT
X
KðkÞ
The summation means an assembly process as same as practice in the FEM, Nedge is the total number of the edges of the whole problem domain O, and KðkÞ is the stiffness matrix associated with edge k given as Z ^ k ÞdO ¼ ðBk ÞT DðB ^ k ÞAk KðkÞ ¼ ðBk ÞT DðB ð31Þ Ok
3.2. Vibration analysis The discretized dynamic equilibrium equation is obtained using the smoothed Galerkin weak form: Z Z T ^ e^ dO ¼ rt duT u€ dO de^ D ð32Þ O
ð26Þ
For boundary edges, the integration domain Ok of edge k is kI kI formed by only one sub-domain, and hence BkI m , Bb and Bs in Eqs. (13), (16) and (22) are obtained using one sub-domain as shown in Fig. 4.
ð30Þ
k¼1
ð24Þ
From Eqs. (12), (17) and (23), we can obtain the smoothed generalized strain in the kth smoothing domain as follows: n oT e^ k ¼ eTmk eTbk eTsk ¼ Bk dk ð25Þ in which h T Bk ¼ Bmk
K¼
O
where r is the density of the material, t is the thickness of the laminated composite plate and u€ the accelerations can be expressed in terms of the nodal displacement di and the shape functions Ni(x) u€ ¼
3 X
Ni ðxÞd€ i
ð33Þ
i¼1
3. Discrete equations
Substituting Eqs. (8), (25) and (33) into Eq. (32) yields KdMd€ ¼ 0
3.1. Bending analysis In the present method, we seek for displacement solution u and the corresponding smoothed strain e^ that satisfy the smoothed Galerkin weak form: Z Z T ^ e^ dO duT f~ dO ¼ 0 de^ D ð27Þ O
O
where f~ represents the external load applied over the problem domain O, and du is the arbitrary variation of the generalized displacement. Substituting Eqs. (8) and (25) into Eq. (27), a set of discretized algebraic system equations can be obtained in the following
ð34Þ
where the smoothed stiffness matrix K has the same expressions as given in Eq. (30), and the mass matrix M is given by MðkÞ ¼ diagfm1 . . .mn g
ð35Þ
where M(k) is sub-matrices of the mass matrix corresponding to edge k, n is the number of the nodes associated with edge k (n ¼3 for boundary edges and n¼4 for inner edges) and mi is the lumped mass at the node i given by n o mi ¼ diag mi mi mi mi t 2 =12 mi t 2 =12 ( rtAk =2 : node i at the end of the edge k mi ¼ ð36Þ 0: node i not at the end of the edge k
822
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
4.2. Bending analysis of laminated plate
A general solution to Eq. (34) can be written as iop t
d ¼ Zp e
ð37Þ
Substituting Eq. (37) in to Eq. (34) yields the eigen equation ðKo2p MÞZp ¼ 0
ð38Þ
where op is the natural frequency associated with the pth mode and Zp is the modal vector.
4. Numerical examples 4.1. Patch test The first numerical example is a homogeneous plate considered for the standard patch test. The mesh the patch with thickness t ¼0.001 is shown in Fig. 5. The plate is subjected to the prescribed displacement and rotations on boundary of the patch (at nodes 1–4) are computed using w ¼ 1þ x þ yþ x2 þxy þ y2 @w ¼ ð1 þ yþ 2xÞ yx ¼ @x @w ¼ ð1 þ x þ 2yÞ yy ¼ @y
ð39Þ
The material properties of patch are E ¼1.0 and n ¼0.25. To satisfy the patch test, the deflection and rotations at any interior nodes computed by numerical method should be in exact agreement the analytic ones given in Eq. (39). To examine the numerical error precisely, an error norm in displacement is defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T exact uPNnode exact u u unum uI unum I I Ed ¼ t I ¼ 1 PNI ð40Þ node exact T uexact I I ¼ 1 uI where the superscript num denotes the displacement vector obtained using numerical methods, exact denotes the exact or analytical solution, and Nnode is the number of total field nodes. Our numerical computation has found the error norm being Ed ¼4.474065656539186 10 15 that is of the order of the machine accuracy. This shows a successful pass of the patch test for present method.
Fig. 5. Mesh for the standard patch test.
4.2.1. Bending of laminated composite plates A cross-ply square plate under uniformly distributed transverse load with different boundary conditions is considered in this subsection, and different length/thickness ratios are used to verify the present method for thick and thin plates. The nine-layer plate symmetrically laminated with [0/90/0/90/0/90/0/90/0] is used. The total thickness of all the 01 and 901 layers is the same, and the layers of the same orientation have the same thickness. The geometry and material properties are given as: L¼20 in, E11 ¼30.0 106 psi, E22 ¼ 7.5 105 psi, G12 ¼4.5 105 psi, G23 ¼3.75 105 psi, G13 ¼4.5 105 psi and u12 ¼u23 ¼u13 ¼ 0.25. The shear correction factors k1 and k2 are taken to be unity. The whole plate is discretized as shown in Fig. 6, and the boundary conditions are given as follows: Simply supported: u ¼ w ¼ yy ¼ 0
at x ¼ 0 and x ¼ L
v ¼ w ¼ yx ¼ 0
at y ¼ 0 and y ¼ L
Clamped: u ¼ v ¼ w ¼ yx ¼ yy ¼ 0
on all the boundaries
ð41Þ
For simply supported plate, the present solutions are compared with 3-node TRIPLT element (15 DOF per node) solutions obtained by Lakshminarayana and Murthy [3]. Numerical results of the maximum transverse displacement at the plate center are given in Fig. 7. The figure indicates the excellent accuracy of these numerical methods. Effect of length/thickness ratios (L/t) on the maximum transverse displacement at the center of the plate is shown in Fig. 8, and no shear locking is observed. Therefore, good results can be obtained using present method regardless of the thickness of the plate. We found also that the present results are very close to the exact solution even if very coarse meshes are used: in Fig. 7, the error is less than 0.3% when only 6 elements per side are used. As a comparison, the corresponding error obtained using TRIPLT element is 2% and more than 12% for extremely thin plate. It is noteworthy that the number of DOF for each node of the present element is 5 that is much smaller than 15 DOFs of the TRIPLT element. The results of DSG method [24] without edge-based smoothing technique are also plotted in Fig. 7. We can find out that the accuracy is observably improved using the edge-based smoothing operation. Fig. 9 shows the numerical results of central deflection for decreasing aspect ratios for a clamped square laminated composite plate. Figs. 8 and 9 show that shear locking can be negligible for both the clamped and simply supported edges.
Fig. 6. Configuration and mesh of the square laminated composite plate.
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
823
L/t = 100
L/t = 10 1.05
1.02 1
1
w/wexact
w/wexact
0.98 0.96 0.94
0.95 Analytical
Analytical
0.92
0.9
Present
Present DSG
DSG
0.9
TRIPLT
TRIPLT
0.85
0.88 4
6
8
10
12
4
Number of elements per side (whole plate)
6
L/t = 200
10
12
L/t = 1000
1.05
1.1 1
w/wexact
1 w/wexact
8
Number of elements per side (whole plate)
0.95
0.9 0.8 Analytical
Analytical
0.9
Present
Present
0.7
DSG
DSG
TRIPLT
TRIPLT
0.85 4
6
8
10
12
Number of elements per side (whole plate)
4
6
8
10
12
Number of elements per side (whole plate)
Fig. 7. Effect of the length/thickness (L/t) on convergence of the central deflection for a simply supported cross-ply square plate subjected to uniformly distributed transverse load.
Fig. 8. Effect of length/thickness (L/t) on the central deflection of a simply supported angle-ply square plate subjected to uniformly distributed transverse load.
Fig. 9. Effect of length/thickness (L/t) on the central deflection of a clamped angle-ply square plate subjected to uniformly distributed transverse load.
4.2.2. Coupled bending and stretching of angle-ply square plate (y/ y) This example considers a two-layered angle-ply laminated graphite/epoxy composite square plate subjected to uniformly distributed load of intensity q. Each layer is of thickness t/2. The plate of length L is simply supported on all four edges. The geometry and material properties of the plate are taken as: L¼ 10 in, t¼0.2 in, E11 ¼40.0 106 psi, E22 ¼1.0 106 psi, G12 ¼5.0 105 psi and u12 ¼0.25. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate is discretized by 16 16 meshes. Numerical
results are obtained for q¼100 psi for 751, 7151, 7251, 7351 and 7451 fiber orientations. For the purpose of testing the element performances, the numerical results obtained from the present method are compared with the 3-node TRIPLT element (15 DOF per node) solutions obtained by Lakshminarayana and Murthy [3] and the hybrid-stress quadrilateral element (MLP3K element) solution obtained by Spilker et al. [25]. Numerical results of maximum values of displacements u, v, w, and bending stress resultants Mxx and Myy are listed in Table 1. From the results, it can be seen that although the present method uses linear triangular elements with
824
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
Table 1 Numerical results of maximum displacements w, u, v and bending stress resultants Mxx, Myy for an angle-ply laminated square plate with two-layer angle-ply laminated graphite/epoxy composites subjected to uniformly distributed load. Quantity
Element
w (at x ¼y ¼L/2)
MLP3K TRIPLT DSG Present
0.609 0.606 0.604 0.608
0.909 0.904 0.900 0.906
0.999 0.992 0.988 0.994
0.960 0.952 0.949 0.954
0.930 0.922 0.920 0.924
u (at x¼ L/2, y ¼0)
MLP3K TRIPLT DSG Present
0.00258 0.00265 0.00256 0.00263
0.00753 0.00773 0.00751 0.00766
0.0112 0.0114 0.0111 0.0113
0.0146 0.0147 0.0143 0.0146
0.0148 0.0148 0.0148 0.0148
v (at x¼ 0, y¼L/2)
MLP3K TRIPLT DSG Present
0.0186 0.0187 0.0185 0.0187
0.0326 0.0327 0.0325 0.0327
0.0249 0.0250 0.0249 0.0249
0.0160 0.0161 0.0164 0.0161
0.0148 0.0148 0.0148 0.0148
Mxx (at x ¼y¼ L/2)
MLP3K TRIPLT DSG Present
Myy (at x¼ y¼L/2)
MLP3K TRIPLT DSG Present
y ¼ 51
1331.0 1327.0 1303.8 1312.4 35.03 34.47 35.48 35.15
only 5 DOF per node, it achieves the same accuracy level as the hybrid-stress quadrilateral element MLP3K and the high order 3-node triangular element TRIPLT with 15 DOF per node. In this example, compared with DSG without edge-based smoothing operation, the present method can obtain more accurate results, especially for bending stresses.
y ¼ 151
1150.0 1150.0 1126.0 1135.1 125.50 124.40 123.21 123.51
y ¼ 351
y ¼451
847.00 851.30 827.27 836.44
569.50 573.10 553.96 560.28
372.20 375.30 364.96 366.04
228.30 228.40 222.52 224.55
307.50 308.80 298.70 301.89
372.20 375.30 364.96 366.04
Table 2 Non-dimensionalized fundamental frequencies o for simply supported four-layer cross-ply [0/90/90/0] square laminated composite plate with different modulus ratio. E11/E22
L/t
Present
DSG [24]
3
5 10
6.526 7.192
6.574 7.249
6.570 7.242
6.360 7.157
10
5 10
8.283 9.810
8.338 9.882
8.271 9.842
8.080 9.670
20
5 10
9.573 12.211
9.633 12.297
9.526 12.218
9.440 12.115
30
5 10
10.341 13.879
10.404 13.974
10.326 13.864
10.238 13.799
40
5 10
10.873 15.142
10.937 15.244
10.854 15.107
10.789 15.068
4.3. Free vibration analysis of laminated plates 4.3.1. Natural frequency analysis of cross-ply laminates A simply supported cross-ply square laminated composite plate is studied with different modulus ratios (E11/E22). The four-layer plate symmetrically laminated with [0/90/90/0] is used and each layer has the same thickness. The following parameters are used in the computation: L¼10.0 m, t¼1 or 2 m, E22 ¼1.0 GPa, G12 ¼G13 ¼0.6E22, G23 ¼0.5E22, u12 ¼u23 ¼u13 ¼0.25 and r ¼1 kg/m3. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate with 12 12 meshes is used in the analysis. Table 2 shows the non-dimensionalized fundamental frequencies, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi o ¼ o11 L2 r=ðE22 t2 Þ, of laminates as functions of modulus ratios (E11/E22). Two kinds of length/thickness ratios are considered. Solutions obtained by Reddy [1] and those of Dai et al. [5] are also provided in Table 2 for comparison purpose. It can be seen that the primary frequency increases with the ratios E11/E22 and L/t. The present results are in a very good agreement with those obtained by other researchers. We can also find that the frequencies obtained using present method is a little smaller than DSG. It indicates that the stiffness of model is softened using the edge-based smoothing operation, which compensates nicely the ‘‘overly-stiff’’ behavior of the FEM model as discussed above. The effect of length/thickness ratios on dimensionless fundapffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mental frequencies (o ¼ o11 L2 r=ðE22 t 2 Þ) is shown in Table 3 with a E11/E22 ratio specified of 40. The present results are compared with solutions given by Reddy and Phan [26], as well as Mallikarjuna and Kant [27]. It is observed that good agreements are attained with Reddy and Phan. For thick plate the present results are in between these two solutions of [26,27]. It is worth pointing out that both Reddy and Mallikarjuna use higher-
y ¼ 251
Reddy [1]
Dai et al. [5]
order shear deformation theory; the present method uses the first-order shear deformation theory for this problem.
4.3.2. Free vibration analysis of angle-ply laminated composite plates Two cases of square angle-ply laminated composite plates are further studied in this section. A four-layer angle-ply [y/ y/y/ y] square plate is first considered. The antisymmetric laminate has the same thickness per ply and is simply supported at all edges. The effect of various fiber orientation angles is investigated. The geometry and material properties are given as: L¼10.0 m, t¼1 m, E11 ¼ 40.0 GPa, E22 ¼1.0 GPa, G12 ¼ G13 ¼0.6E22, G23 ¼0.5E22, u12 ¼u23 ¼u13 ¼0.25 and r ¼1 kg/m3. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate with 16 16 meshes is used in the analysis. The exact solutions given by Bert and Chen [28], the finite strip solutions given by Wang and Lin [29], and other finite element solutions given by Mallikarjuna and Kant [27], Reddy [2], are used for comparison. The fundamental frequencies normalized by o ¼ o11 L2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r=ðE22 t2 Þ versus orientation angle y are tabulated in Table 4. It
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
825
Table 3 Effect of aspect ratio L/t on the non-dimensionalized fundamental frequencies o for a simply supported four-layer cross-ply [0/90/90/0] square laminated composite plate. L/t
4
5
10
20
50
100
Present DSG [24] Reddy and Phan [26] Mallikarjuna and Kant [27]
9.406 9.459 9.497 9.258
10.873 10.937 10.989 10.740
15.142 15.244 15.270 15.090
17.605 17.733 17.668 17.637
18.588 18.727 18.606 18.669
18.746 18.885 18.755 18.835
Table 4 Non-dimensional fundamental frequencies o for a simply supported four-layer angle-ply [y/ y/y/ y] square laminated plate with varying fiber orientations.
y
51
151
301
451
Present DSG [24] Bert and Chen [28] Wang and Lin [29] Mallikarjuna and Kant [27] Reddy [2]
14.686 14.743 14.740 14.853 – –
15.540 15.604 15.610 15.759 – –
17.544 17.622 17.630 17.702 17.513 –
18.323 18.405 18.460 18.263 18.320 18.259
Table 5 Non-dimensional fundamental frequencies on for a simply supported ten-layer angle-ply [y/ y/y] square laminated plate with varying aspect L/t. L/t
100 20 10 6.67 5 4 3.33
y ¼ 151
y ¼ 301
domains. The rectangle rule is employed for performing the numerical integration over each smoothing domain. Based on the results from extensive numerical evaluations, the following features of the present method are found. (1) Present ES-FEM uses 3-node triangular element with five degrees of freedom per node. There is no extra sampling point introduced to evaluate the stiffness matrix in present formulation. The formulation is very simple, and its performance is significantly improved by using an edge-based smoothing operation. (2) The solutions obtained using 3-noded triangular elements in the present study are very stable and accurate in comparison to other complex plate elements. (3) Satisfactory accuracy and rate of convergence can be obtained for both thin and thick the laminated composite plates. (4) The present element can be easily extended to quadrangular element and general n-sided (n44) polygonal elements.
y ¼451
Present
Noor and Burton [4]
Present
Noor and Burton [4]
Present
Noor and Burton [4]
0.001323 0.03193 0.1160 0.2303 0.3586 0.4927 0.6290
0.001328 0.03201 0.1162 0.2304 0.3588 0.4934 0.6307
0.001503 0.03611 0.1299 0.2548 0.3925 0.5341 0.6762
0.001510 0.03620 0.1296 0.2532 0.3889 0.5286 0.6692
0.001586 0.03796 0.1353 0.2630 0.4023 0.5445 0.6866
0.001595 0.03808 0.1351 0.2617 0.3993 0.5400 0.6810
can be found that the present results agree well with those compared solutions. The maximum deviation of the present solutions from the exact solutions, which are given by Bert and Chen [28], is only 0.72%. A ten-layer angle-ply [y/ y/]5 laminated composite square plate of length L¼1.0 m is used for studying the effect of length/ thickness ratios. The laminate has the same thickness per layer and is simply supported at all edges. Seven different thickness, t ¼0.01, 0.05, 0.10, 0.15, 0.20, 0.25 and 0.30 m, are considered in this study. The material properties are given as: E11 ¼15.0 GPa, E22 ¼1.0 GPa, G12 ¼G13 ¼ 0.5E22, G23 ¼0.35E22, u12 ¼ u13 ¼0.30, u23 ¼0.49 and r ¼1 kg/m3. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate with 20 20 meshes is used in the analysis. The elasticity solutions obtained by Noor and Burton [4] are used as p the reference solutions. The solutions ffiffiffiffiffiffiffiffiffiffiffiffi ffi normalized by o ¼ o11 t r=E22 are shown in Table 5. It can be found that the present solutions agree well with the reference solutions for all the length/thickness ratios tested here.
5. Conclusion Bending and vibration analysis of laminated composite plates has been studied based on an edge-based smoothing technique for using simple 3-node triangular elements. In the present formulation, the smoothed Galerkin weak form is adopted to build the discretized system equations in edge-based smoothing
Acknowledgments The support of National Science Foundation of China (11002053), National Basic Research Program of China (2010CB328005), Program for Changjiang Scholar and Innovative Research Team in University, and Science Fund of State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body (51075005) are gratefully acknowledged.
References [1] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. CRC press; 2004. [2] Reddy JN. Static free vibration of antisymmetric, angle-ply laminated plates including transverse shear deformation by the finite element method. Journal of Sound and Vibration 1979;66(4):565–76. [3] Lakshminarayana HV, Murthy SS. A shear-flexible triangular finite element model for laminated composite plates. International Journal for Numerical Methods in Engineering 1984;20(4):591–623. [4] Noor AK, Burton WS. Assessment of computational models for multilayered anisotropic plates. Composite Structures 1990;14(3):233–65. [5] Dai KY, Liu GR, Lim KM, Chen XL. A mesh-free method for static and free vibration analysis of shear deformable laminated composite plates. Journal of Sound and Vibration 2004;269(3–5):633–52. [6] Setoodeh AR, Karami G. Static free vibration and buckling analysis of anisotropic thick laminated composite plates on distributed and point elastic supports using a 3-D layer-wise FEM. Engineering Structures 2004;26:211–20. [7] Lee WH, Han SC. Free and forced vibration analysis of laminated composite plates and shells using 9-node assumed strain shell element. Computational Mechanics 2006;39(1):41–58. [8] Liu GR. A generalized gradient smoothing technique and smoothed bilinear form for Galerkin formulation of a wide class of computational methods. International Journal of Computational Methods 2008;5(2):199–236. [9] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin meshfree methods. International Journal for Numerical Methods in Engineering 2001;50:435–66. [10] Liu GR, Zhang GY, Dai KY, Wang YY, Zhong ZH, Li GY, et al. A linearly conforming point interpolation method (LC-PIM) for 2D solid mechanics problems. International Journal of Computational Methods 2005;2:645–65. [11] Liu GR, Li Y, Dai KY, Luan MT, Xue WA. Linearly conforming radial point interpolation method for solid mechanics problems. International Journal of Computational Methods 2006;3(4):401–28. [12] Liu GR, Zhang GY. Upper bound solution to elasticity problems: a unique property of the linearly conforming point interpolation method (LC-PIM). International Journal for Numerical Methods in Engineering 2008;74(7):1128–61.
826
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
[13] Liu GR, Zhao X, Dai KY, Zhong ZH, Li GY, Han X. Static and free vibration analysis of laminated composite plates using the conforming radial point interpolation method. Composites Science and Technology 2008;68(2):354–66. [14] Puso MA, Solberg J. A stabilized nodally integrated tetrahedral. International Journal for Numerical Methods in Engineering 2006;67:841–67. [15] Puso MA, Chen JS, Zywicz E, Han X, Elmer W. Meshfree and finite element nodal integration methods. International Journal for Numerical Methods in Engineering 2008;74:416–46. [16] Nagashima T. Node-by-node meshless approach and its applications to structural analyses. International Journal for Numerical Methods in Engineering 1999;46:341–85. [17] Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. Computational Mechanics 2007;39(6):859–77. [18] Liu GR, Nguyen TT, Dai KY, Lam KY. Theoretical aspects of the smoothed finite element method (SFEM). International Journal for Numerical Methods in Engineering 2007;71(8):902–30. [19] Cui XY, Liu GR, Li GY, Zhao X, Nguyen TT, Sun GY. A smoothed finite element method (SFEM) for linear and geometrically nonlinear analysis of plates and shells. CMES: Computer Modeling in Engineering & Science 2008;28(2):109–26. [20] Nguyen-Xuan H, Rabczuk T, Bordas S, Debongnie JF. A smoothed finite element method for plate analysis. Computer Methods in Applied Mechanics and Engineering 2008;197:1184–203. [21] Nguyen-Xuan H, Nguyen-Thoi T. A stabilized smoothed finite element method for free vibration analysis of Mindlin–Reissner plates. Communications in Numerical Methods in Engineering 2009;25:882–906.
[22] Dai KY, Liu GR, Nguyen. TT. An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics. Finite Elements in Analysis and Design 2007;43:847–60. [23] Liu GR, Nguyen TT, Lam. KY. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analysis. Journal of Sound and Vibration 2009;320:1100–30. [24] Bletzinger KU, Bischoff M, Ramm. E. A unified approach for shear-locking-free triangular and rectangular shell finite elements. Computers & Structures 2000;75(3):321–34. [25] Spilker RL, Chou SC, Orringer O. Alternate hybrid-stress elements for analysis of multilayer composite plates. Journal of Composite Materials 1977;11(1):51–70. [26] Reddy JN, Phan ND. Stability and vibration of isotropic, orthotropic, and laminated plates using a higher-order shear deformation theory. Journal of Sound and Vibration 1985;98(2):157–70. [27] Mallikarjuna T Kant. Free vibration of symmetrically laminated plates using a higher order theory with finite element technique. International Journal for Numerical Methods in Engineering 1989;28(8):1875–89. [28] Bert CW, Chen TLC. Effect of shear deformation on vibration of antisymmetric angle-ply laminated rectangular plates. International Journal of Solids and Structures 1978;14(6):465–73. [29] Wang WJ, Lin KC. Free vibration of laminated plates using a finite strip method based on a higher-order plate theory. Computers & Structures 1994;53(6): 1281–9.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Bending and vibration responses of laminated composite plates using an edge-based smoothing technique X.Y. Cui a, G.R. Liu b, G.Y. Li a,n a b
State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, PR China School of Aerospace Systems, University of Cincinnati, Cincinnati, OH 45221-0070, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 November 2009 Accepted 16 December 2010 Available online 25 February 2011
In this paper, bending and vibration analysis of laminated composite plates is carried out using a novel triangular composite plate element based on an edge-based smoothing technique. The present formulation is based on the first-order shear deformation theory, and the discrete shear gap (DSG) method is employed to mitigate the shear locking. The smoothed Galerkin weak form is adopted to obtain the discretized system equations, and edge-based smoothing domains are used for the numerical integration to improve the accuracy and the convergence rate of the method. The present formulation is coded and used to solve various example problems of bending and free vibration of laminated composite plates. It is found that the present method can provide excellent results with a wide range of thickness and is free of shear locking. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Numerical methods Vibration FEM Techniques-FEM Laminated composite plates
1. Introduction The bending and vibration analysis of laminated composite plates has become more and more important in the past several decades, due to rapidly increased application of composite materials in structure systems in civil, mechanical and aerospace engineering. The search for simple, efficient and inexpensive numerical methods applicable to both thin and thick laminated composite plates has been receiving much interest. Various laminated plate theories have been introduced for analysis of thin to thick laminated plates such as the classical laminated plate theory (CLPT), the first-order shear deformation theory (FSDT) and the higher-order shear deformation theory (HSDT) [1]. Using these theories, some existing numerical methods have been adopted to compute bending deformation and natural frequencies of composite laminates. Reddy [2] developed the free vibration analysis of antisymmetric, angle-ply laminated plates including transverse shear deformation using the finite element method. A 3-node higher-order triangular plate element TRIPLT was formulated by Lakshminarayana and Murthy [3]. This element has 15 degrees of freedom (DOF) containing displacements and rotations along with their first derivatives per node. Noor and Burton [4] gave an assessment of computational models for multilayered anisotropic plates. Dai et al. [5] studied the static and free vibration analysis based on HSDT using element free Galerkin
n
Corresponding author. Tel.: +86 731 8821717; fax: + 86 731 8822051. E-mail address: [email protected] (G.Y. Li).
0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.01.007
(EFG) method. Setoodeh and Karami [6] investigated the characteristics of anisotropic thick laminated composite plates using a 3-D layer-wise FEM. Lee and Han [7] developed a 9-node assumed strain element for free and forced vibration analysis of laminated composite plates and shells. Although large numbers of laminated plate elements have been developed, most of them have a complex form. Displacement-based triangular elements have a very simple form as they only consider the displacement and rotation degrees of freedom at nodes and without additional parameters. However, these elements generally have a low accuracy because of the inherent characteristic of overstiff phenomenon from the standard weak formulation. In order to ‘‘soften’’ the system and allow the use of discontinuous shape functions, Liu [8] proposed a general formulation of smoothed Galerkin weak form using the strain smoothing technique [9]. Efficient numerical methods have been developed based on both FEM and meshfree settings using different smoothing domains that can be node-based cell-based can edge-based. Using the point interpolation method and the node-based smoothing operations, a node-based smoothed point interpolation method (NS-PIM) were developed using polynomial basis functions (NS-PIM [10]) and radial basis functions (LC-RPIM [11]). Liu and Zhang [12] proved that the NS-PIM is variationally consistent, and can provide much better stress results. More importantly it can provide upper bound solution in energy norm. Liu et al. [13] has found that the LC-RPIM needs a high order interpolation to obtain satisfactory solutions for laminated composite plates. However, the node-based smoothed models are ‘‘overly-soft’’, which leads to temporal instability [14–16]. The instability can be clearly shown as spurious non-zero energy modes
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
in free vibration analyses. Using the FEM shape functions and strain smoothing operations based on cells that are obtained by dividing the quadrilateral elements (Q4), a smoothed finite element method (SFEM) was first formulated by Liu and co-workers [17]. Liu et al. [18] gave detailed theoretical aspects including stability, bound property and convergence about SFEM, and Cui et al. [19] extended the SFEM for linear and geometrically nonlinear analysis of plates and shells. The SFEM idea was also applied to the MITC4 element by Nguyen-Xuan et al. [20] for linear and static analysis and free vibration analysis [21] of plates. The formulation of SFEM is found much flexible and robust than the FEM, and was further extended to general n-sided polygonal elements [22]. However, the SFEM does not give any advantage, when triangular elements (T3) are used, because the SFEM-T3 is in fact exactly the same as the FEM-T3. Recently, an edge-based smoothed finite element (ES-FEM) is proposed by constructing edge-based smoothing domain using the smoothed Galerkin weak form. The ES-FEM exhibit superconvergence properties in both energy norm and displacement norm and can produce ultra-accurate solutions in 2D solid mechanics problems [23]. In this work, the ES-FEM is further extended to solve bending and vibration problems of laminated composite plates using simple 3-node triangular elements that can be generated automatically. The present formulation is based on the first-order shear deformation theory, and shear locking effect is suppressed using the discrete shear gap method [24]. In the present method, the smoothing domains associated with the edges of the triangles are formed based on the triangular element mesh. The compatible strains in the smoothing domains are smoothed using the simple rectangle rule of integration, and the system stiffness matrix is obtained by numerical integration based on smoothing domains. The present method has a very simple form and without using any internal parameters in the process of building the stiffness matrix. Using the present ES-FEM method, excellent results have been obtained for both bending and vibration problems of laminated composite plates.
819
of the transverse normals, to the undeformed mid-plane, in the xz and yz planes, respectively. The strain vector e can be written in terms of the mid-plane deformations: ( ) 0 em zeb e ¼ fexx , eyy , gxy , gxz , gyz gT ¼ þ þ ð2Þ e 0 0 s where em is the membrane strain, eb is the bending strain and es is the shear strain, which can be given by 9 9 8 8 @y @u0 x 8 9 > > > > > > > > @x @x > > > > 0 = = < < < @w = @yy @v0 @x yx , eb ¼ , es ¼ @w0 ð3Þ em ¼ @y @y > > > > : @y yy ; > > > > > > ; ; : @u0 þ @v0 > : @yx þ @yy > @y
@x
@y
@x
The constitutive relationship of laminated composite plates can be expressed as
r^ ¼ D^ e^
ð4Þ
where 2
r^ ¼ f r^ Tm r^ Tb r^ Ts gT , e^ ¼ f eTm eTb e
T T s g ,
3
A
B
0
^ ¼6 D 4B 0
D
7 05
0
S
ð5Þ
^ m ¼ fNx ,Ny ,Nxy gT is the membrane force vector in the midwhere r ^ b ¼ fMx ,My ,Mxy gT is the bending moment vector and plane, r r^ s ¼ fQx , Qy gT is the transverse shear force vector. The membrane stiffness (A), the bending stiffness (D), the bending–extensional coupling stiffness (B) and the transverse shear stiffness (S) are defined as 2 3 2 3 A11 A12 A16 B11 B12 B16 6A 7 6 7 A ¼ 4 12 A22 A26 5, B ¼ 4 B12 B22 B26 5 A16 A26 A66 B16 B26 B66 2 3 " # D11 D12 D16 pffiffiffiffiffiffiffiffiffiffiffi k1 S44 k1 k2 S45 6D 7 D D p ffiffiffiffiffiffiffiffiffiffiffi ð6Þ S¼ D ¼ 4 12 22 26 5, k1 k2 S45 k2 S55 D16 D26 D66
2. Basic formulations in which k1 and k2 are shear correction factors. The stiffness in Eq. (6) can be defined as
2.1. Geometry and kinematics of laminated composite plate A typical laminated composite plate with thickness t and Nl layers is shown in Fig. 1. Based on the first-order shear deformation laminated plate theory [1], the displacements field can be expressed by uðx,y,zÞ ¼ u0 ðx,yÞ þzyx ðx,yÞ vðx,y,zÞ ¼ v0 ðx,yÞ þ zyy ðx,yÞ wðx,y,zÞ ¼ w0 ðx,yÞ
ð1Þ
in which u0, v0 and w0 are the displacements of the mid-plane of the plate in the x, y and z directions, yx and yy denote the rotations
Aij ¼ Dij ¼
Nl X
ðkÞ
Q ij ðtk tk þ 1 Þ,
k¼1 Nl X
1 ðkÞ Q ðt 3 t 3 Þ, 3 k ¼ 1 ij k k þ 1
Bij ¼
Nl 1X ðkÞ Q ðt 2 t 2 Þ 2 k ¼ 1 ij k k þ 1
Sij ¼
Nl X
ðkÞ
Q ij ðtk tk þ 1 Þ
ð7Þ
k¼1
where Aij, Bij and Dij are defined for i, j ¼1, 2, 6, whereas Sij is defined for i, j ¼4, 5. tk and tk + 1 denote the distances from the plate reference mid-plane to the outer and inner surfaces of the kth layer, respectively, as shown in Fig. 2. Nl is the total number of ðkÞ
layers in the laminated plate, and Q ij (i, j¼1, 2, 4, 5, 6) are related to the engineering constants for the kth layer and can refer to the work by Reddy [1].
t
2.2. Edge-based strain smoothing technique
Fig. 1. A typical laminated composite plate with thickness t and Nl layers and its coordinate system.
In the present method, the domain discretization is based on triangular elements which are same as those used in standard FEM. The problem domain O is divided into Ne triangular elements with a total of Nedge edges. The smoothing domain for each edge is formed by sequentially connecting two end points of the edge and centroids of its surrounding triangles to perform the smoothing operation, such that O ¼ O1 [ O2 [ [ ONedge and Oi\Oj ¼| (iaj, i¼1, y, Nedge, j¼ 1, y, Nedge). A strain smoothing
820
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
Fig. 3. The integration domain Ok of edge k is assembled by two sub-domains, Ok1 and Ok2. The sub-domain Ok1 is from the element e1, and the sub-domain Ok2 is form the element e2, xCk1 and xCk2 are the geometric centers of the sub-domains. Fig. 2. Cross-sectional view of the laminated plate.
technique is used over each smoothing domains and the integration is performed based on the smoothing domains associated with the edges, instead of based on the elements as in the standard FEM. The displacement at any point in the problem domain is interpolated using the element shape functions, which are created as same as in the linear standard FEM. The generalized mid-plane displacements vector u¼{u0, v0, w0, yx, yy}T consisting of both displacements and rotations can be then approximated as uðxÞ ¼
3 X
Ni ðxÞdi
ð8Þ
i¼1
where di ¼ {u0i, v0i, wi, yxi, yyi}T is the displacement of the node i, and Ni is the diagonal matrix of shape functions n
zfflfflfflfflfflfflffl}|fflfflfflfflfflfflffl{ Ni ¼ diagfg fNi ,. . .,Ni g
ð9Þ
in which n ¼5 is the degree of freedom in one node and Ni is the linear shape function. Introducing the strain smoothing operation [8], the strain in the smoothing domain Ok is given by Z 1 ek ¼ e^ ðxÞdO ð10Þ Ak Ok where e^ ðxÞ is the generalized compatible strain given in Eq. (4) and Ak is the area of the smoothing domain Ok. Since linear shape functions are used in the present method, the compatible strain e^ ðxÞ is a constant in each smoothing subdomain OkI. Therefore, the integration in Eq. (10) can be performed exactly using the simple rectangle rule of integration:
ek ¼
1 X AkI e^ ðxCkI Þ Ak
ð11Þ
where xCkI is the centroid of the smoothing sub-domain OkI, and AkI is the area of the smoothing sub-domain OkI. It is clear that all one needs to do is to sample the compatible strains at the geometric center of the smoothing sub-domain, which can be done very efficiently. As shown in Fig. 3, for interior edges, the integration domain Ok of edge k is formed by two sub-domains, Ok1 and Ok2. The subdomains Ok1 and Ok2 are from the element e1 and the element e2, respectively. The nodes of elements e1 and e2 are i, j, m and j, n, m, respectively. Substituting Eq. (8) into Eq. (3) and using the strain smoothing technique, the smoothed membrane strain can be rewritten as 1 emk ¼ ½Ak1 em ðxCk1 Þ þ Ak2 em ðxCk1 Þ Ak A A ðx Þd þ k2 Bk2 ðx Þd ¼ Bmk dk ¼ k1 Bk1 Ak m Ck1 k1 Ak m Ck2 k2
where dk1 and dk2 are the nodal displacement vectors of the element e1 and the element e2, respectively, dk are the displacements of the nodes associated with edge k, and BkI m ðxCkI Þ are the compatible strain matrices of the smoothing sub-domain OkI (I ¼1, 2), which are given as h i kI BkI BkI ðI ¼ 1,2Þ BkI m2 m3 , m ðxCk1 Þ ¼ Bm1 2 @N ðx Þ 3 i CkI 0 0 0 0 @x 6 7 @Ni ðxCkI Þ 6 0 0 0 07 BkI ð13Þ 7, ði ¼ 1,2,3Þ @y mi ¼ 6 4 5 @Ni ðxCkI Þ @Ni ðxCkI Þ 0 0 0 @y @x From Eq. (12), the smoothed membrane strain matrix Bmk can be written as Bmk ¼
Ak1 k1 Ak2 k2 B þ B Ak m Ak m
ð14Þ
Note that the sign ‘+ ’ denotes assembly but not summation here. Similarly, we can obtain the bending strain matrix Bbk as follows: Bbk ¼
Ak1 k1 Ak2 k2 B þ B Ak b Ak b
ð15Þ
in which h kI BkI b ðxCk1 Þ ¼ Bb1 2 0 0 0 6 60 0 0 kI Bbi ¼ 6 4 0 0 0
BkI b2
i BkI b3 ,
@Ni ðxCkI Þ @x
0 @Ni ðxCkI Þ @y
0
ðI ¼ 1,2Þ 3 7
@Ni ðxCkI Þ 7 7, @y @Ni ðxCkI Þ @x
ði ¼ 1,2,3Þ
5
The smoothed bending strain can be given as follows:
ebk ¼ Bbk dk
ð17Þ
In order to eliminate the shear locking, the ‘‘discrete shear gap’’ (DSG) method [24] is used. In each triangular element, the shear strain here can be computed by
gxz ¼
3 X @Ni i¼1
gyz ¼
@x
3 X @Ni i¼1
@y
3 X @Ni
Dwxi þ
i¼1
@x
3 X @Ni
Dwxi þ
i¼1
@y
Dwyi Dwyi
ð18Þ
where Dwxi and Dwyi are the discrete shear gaps at the node i, they are given by
Dwx1 ¼ Dwx3 ¼ Dwy1 ¼ Dwy2 ¼ 0 ð12Þ
ð16Þ
1 2
1 2
Dwx2 ¼ ðw2 w1 Þ þ aðyx1 þ yx2 Þ þ bðyy1 þ yy2 Þ
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
1 2
1 2
Dwy3 ¼ ðw3 w1 Þ þ cðyx1 þ yx3 Þ þ dðyy1 þ yy3 Þ
821
ð19Þ
in which a ¼ x2 x1 ,
b ¼ y2 y1
c ¼ x3 x1 ,
d ¼ y3 y1
ð20Þ
From Eqs. (18)–(20), the shear strain es in each element can be written as n oT es ¼ gxz gyz ¼ Bs de ð21Þ where de is the nodal displacement of the triangular element, Bs is the shear strain matrix given as Bs ¼ Bs1 Bs2 Bs3 " # 1 0 0 bd Ae 0 Bs1 ¼ 2Ae 0 0 ca 0 Ae " # ad bd 1 0 0 d 2 2 Bs2 ¼ 2Ae 0 0 c ac bc 2 2 " # bc bd 0 0 b 1 2 2 Bs3 ¼ ð22Þ ac ad 2Ae 0 0 a 2 2
Fig. 4. The integration domain Ok of edge k is assembled by only one domain.
matrix form: Kdf ¼ 0
ð28Þ
where f is the force vector defined as Z f ¼ NðxÞf~ dO
ð29Þ
O
and the smoothed stiffness matrix K are expressed as Nedge
For the smoothing domain Ok assembled by sub-domains Ok1 and Ok2, The shear strain matrices of the correlative triangular k2 element Bk1 s and Bs are computed from Eq. (22). The smoothed k shear strain es of the domain Ok can be written as
esk ¼ Bsk dk
ð23Þ
where Bsk is assembled by Bsk ¼
Ak1 k1 Ak2 k2 B þ B Ak s Ak s
Bk1 s
and
Bk2 s
as follows:
T
Bbk
T
Bsk
iT
X
KðkÞ
The summation means an assembly process as same as practice in the FEM, Nedge is the total number of the edges of the whole problem domain O, and KðkÞ is the stiffness matrix associated with edge k given as Z ^ k ÞdO ¼ ðBk ÞT DðB ^ k ÞAk KðkÞ ¼ ðBk ÞT DðB ð31Þ Ok
3.2. Vibration analysis The discretized dynamic equilibrium equation is obtained using the smoothed Galerkin weak form: Z Z T ^ e^ dO ¼ rt duT u€ dO de^ D ð32Þ O
ð26Þ
For boundary edges, the integration domain Ok of edge k is kI kI formed by only one sub-domain, and hence BkI m , Bb and Bs in Eqs. (13), (16) and (22) are obtained using one sub-domain as shown in Fig. 4.
ð30Þ
k¼1
ð24Þ
From Eqs. (12), (17) and (23), we can obtain the smoothed generalized strain in the kth smoothing domain as follows: n oT e^ k ¼ eTmk eTbk eTsk ¼ Bk dk ð25Þ in which h T Bk ¼ Bmk
K¼
O
where r is the density of the material, t is the thickness of the laminated composite plate and u€ the accelerations can be expressed in terms of the nodal displacement di and the shape functions Ni(x) u€ ¼
3 X
Ni ðxÞd€ i
ð33Þ
i¼1
3. Discrete equations
Substituting Eqs. (8), (25) and (33) into Eq. (32) yields KdMd€ ¼ 0
3.1. Bending analysis In the present method, we seek for displacement solution u and the corresponding smoothed strain e^ that satisfy the smoothed Galerkin weak form: Z Z T ^ e^ dO duT f~ dO ¼ 0 de^ D ð27Þ O
O
where f~ represents the external load applied over the problem domain O, and du is the arbitrary variation of the generalized displacement. Substituting Eqs. (8) and (25) into Eq. (27), a set of discretized algebraic system equations can be obtained in the following
ð34Þ
where the smoothed stiffness matrix K has the same expressions as given in Eq. (30), and the mass matrix M is given by MðkÞ ¼ diagfm1 . . .mn g
ð35Þ
where M(k) is sub-matrices of the mass matrix corresponding to edge k, n is the number of the nodes associated with edge k (n ¼3 for boundary edges and n¼4 for inner edges) and mi is the lumped mass at the node i given by n o mi ¼ diag mi mi mi mi t 2 =12 mi t 2 =12 ( rtAk =2 : node i at the end of the edge k mi ¼ ð36Þ 0: node i not at the end of the edge k
822
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
4.2. Bending analysis of laminated plate
A general solution to Eq. (34) can be written as iop t
d ¼ Zp e
ð37Þ
Substituting Eq. (37) in to Eq. (34) yields the eigen equation ðKo2p MÞZp ¼ 0
ð38Þ
where op is the natural frequency associated with the pth mode and Zp is the modal vector.
4. Numerical examples 4.1. Patch test The first numerical example is a homogeneous plate considered for the standard patch test. The mesh the patch with thickness t ¼0.001 is shown in Fig. 5. The plate is subjected to the prescribed displacement and rotations on boundary of the patch (at nodes 1–4) are computed using w ¼ 1þ x þ yþ x2 þxy þ y2 @w ¼ ð1 þ yþ 2xÞ yx ¼ @x @w ¼ ð1 þ x þ 2yÞ yy ¼ @y
ð39Þ
The material properties of patch are E ¼1.0 and n ¼0.25. To satisfy the patch test, the deflection and rotations at any interior nodes computed by numerical method should be in exact agreement the analytic ones given in Eq. (39). To examine the numerical error precisely, an error norm in displacement is defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T exact uPNnode exact u u unum uI unum I I Ed ¼ t I ¼ 1 PNI ð40Þ node exact T uexact I I ¼ 1 uI where the superscript num denotes the displacement vector obtained using numerical methods, exact denotes the exact or analytical solution, and Nnode is the number of total field nodes. Our numerical computation has found the error norm being Ed ¼4.474065656539186 10 15 that is of the order of the machine accuracy. This shows a successful pass of the patch test for present method.
Fig. 5. Mesh for the standard patch test.
4.2.1. Bending of laminated composite plates A cross-ply square plate under uniformly distributed transverse load with different boundary conditions is considered in this subsection, and different length/thickness ratios are used to verify the present method for thick and thin plates. The nine-layer plate symmetrically laminated with [0/90/0/90/0/90/0/90/0] is used. The total thickness of all the 01 and 901 layers is the same, and the layers of the same orientation have the same thickness. The geometry and material properties are given as: L¼20 in, E11 ¼30.0 106 psi, E22 ¼ 7.5 105 psi, G12 ¼4.5 105 psi, G23 ¼3.75 105 psi, G13 ¼4.5 105 psi and u12 ¼u23 ¼u13 ¼ 0.25. The shear correction factors k1 and k2 are taken to be unity. The whole plate is discretized as shown in Fig. 6, and the boundary conditions are given as follows: Simply supported: u ¼ w ¼ yy ¼ 0
at x ¼ 0 and x ¼ L
v ¼ w ¼ yx ¼ 0
at y ¼ 0 and y ¼ L
Clamped: u ¼ v ¼ w ¼ yx ¼ yy ¼ 0
on all the boundaries
ð41Þ
For simply supported plate, the present solutions are compared with 3-node TRIPLT element (15 DOF per node) solutions obtained by Lakshminarayana and Murthy [3]. Numerical results of the maximum transverse displacement at the plate center are given in Fig. 7. The figure indicates the excellent accuracy of these numerical methods. Effect of length/thickness ratios (L/t) on the maximum transverse displacement at the center of the plate is shown in Fig. 8, and no shear locking is observed. Therefore, good results can be obtained using present method regardless of the thickness of the plate. We found also that the present results are very close to the exact solution even if very coarse meshes are used: in Fig. 7, the error is less than 0.3% when only 6 elements per side are used. As a comparison, the corresponding error obtained using TRIPLT element is 2% and more than 12% for extremely thin plate. It is noteworthy that the number of DOF for each node of the present element is 5 that is much smaller than 15 DOFs of the TRIPLT element. The results of DSG method [24] without edge-based smoothing technique are also plotted in Fig. 7. We can find out that the accuracy is observably improved using the edge-based smoothing operation. Fig. 9 shows the numerical results of central deflection for decreasing aspect ratios for a clamped square laminated composite plate. Figs. 8 and 9 show that shear locking can be negligible for both the clamped and simply supported edges.
Fig. 6. Configuration and mesh of the square laminated composite plate.
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
823
L/t = 100
L/t = 10 1.05
1.02 1
1
w/wexact
w/wexact
0.98 0.96 0.94
0.95 Analytical
Analytical
0.92
0.9
Present
Present DSG
DSG
0.9
TRIPLT
TRIPLT
0.85
0.88 4
6
8
10
12
4
Number of elements per side (whole plate)
6
L/t = 200
10
12
L/t = 1000
1.05
1.1 1
w/wexact
1 w/wexact
8
Number of elements per side (whole plate)
0.95
0.9 0.8 Analytical
Analytical
0.9
Present
Present
0.7
DSG
DSG
TRIPLT
TRIPLT
0.85 4
6
8
10
12
Number of elements per side (whole plate)
4
6
8
10
12
Number of elements per side (whole plate)
Fig. 7. Effect of the length/thickness (L/t) on convergence of the central deflection for a simply supported cross-ply square plate subjected to uniformly distributed transverse load.
Fig. 8. Effect of length/thickness (L/t) on the central deflection of a simply supported angle-ply square plate subjected to uniformly distributed transverse load.
Fig. 9. Effect of length/thickness (L/t) on the central deflection of a clamped angle-ply square plate subjected to uniformly distributed transverse load.
4.2.2. Coupled bending and stretching of angle-ply square plate (y/ y) This example considers a two-layered angle-ply laminated graphite/epoxy composite square plate subjected to uniformly distributed load of intensity q. Each layer is of thickness t/2. The plate of length L is simply supported on all four edges. The geometry and material properties of the plate are taken as: L¼ 10 in, t¼0.2 in, E11 ¼40.0 106 psi, E22 ¼1.0 106 psi, G12 ¼5.0 105 psi and u12 ¼0.25. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate is discretized by 16 16 meshes. Numerical
results are obtained for q¼100 psi for 751, 7151, 7251, 7351 and 7451 fiber orientations. For the purpose of testing the element performances, the numerical results obtained from the present method are compared with the 3-node TRIPLT element (15 DOF per node) solutions obtained by Lakshminarayana and Murthy [3] and the hybrid-stress quadrilateral element (MLP3K element) solution obtained by Spilker et al. [25]. Numerical results of maximum values of displacements u, v, w, and bending stress resultants Mxx and Myy are listed in Table 1. From the results, it can be seen that although the present method uses linear triangular elements with
824
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
Table 1 Numerical results of maximum displacements w, u, v and bending stress resultants Mxx, Myy for an angle-ply laminated square plate with two-layer angle-ply laminated graphite/epoxy composites subjected to uniformly distributed load. Quantity
Element
w (at x ¼y ¼L/2)
MLP3K TRIPLT DSG Present
0.609 0.606 0.604 0.608
0.909 0.904 0.900 0.906
0.999 0.992 0.988 0.994
0.960 0.952 0.949 0.954
0.930 0.922 0.920 0.924
u (at x¼ L/2, y ¼0)
MLP3K TRIPLT DSG Present
0.00258 0.00265 0.00256 0.00263
0.00753 0.00773 0.00751 0.00766
0.0112 0.0114 0.0111 0.0113
0.0146 0.0147 0.0143 0.0146
0.0148 0.0148 0.0148 0.0148
v (at x¼ 0, y¼L/2)
MLP3K TRIPLT DSG Present
0.0186 0.0187 0.0185 0.0187
0.0326 0.0327 0.0325 0.0327
0.0249 0.0250 0.0249 0.0249
0.0160 0.0161 0.0164 0.0161
0.0148 0.0148 0.0148 0.0148
Mxx (at x ¼y¼ L/2)
MLP3K TRIPLT DSG Present
Myy (at x¼ y¼L/2)
MLP3K TRIPLT DSG Present
y ¼ 51
1331.0 1327.0 1303.8 1312.4 35.03 34.47 35.48 35.15
only 5 DOF per node, it achieves the same accuracy level as the hybrid-stress quadrilateral element MLP3K and the high order 3-node triangular element TRIPLT with 15 DOF per node. In this example, compared with DSG without edge-based smoothing operation, the present method can obtain more accurate results, especially for bending stresses.
y ¼ 151
1150.0 1150.0 1126.0 1135.1 125.50 124.40 123.21 123.51
y ¼ 351
y ¼451
847.00 851.30 827.27 836.44
569.50 573.10 553.96 560.28
372.20 375.30 364.96 366.04
228.30 228.40 222.52 224.55
307.50 308.80 298.70 301.89
372.20 375.30 364.96 366.04
Table 2 Non-dimensionalized fundamental frequencies o for simply supported four-layer cross-ply [0/90/90/0] square laminated composite plate with different modulus ratio. E11/E22
L/t
Present
DSG [24]
3
5 10
6.526 7.192
6.574 7.249
6.570 7.242
6.360 7.157
10
5 10
8.283 9.810
8.338 9.882
8.271 9.842
8.080 9.670
20
5 10
9.573 12.211
9.633 12.297
9.526 12.218
9.440 12.115
30
5 10
10.341 13.879
10.404 13.974
10.326 13.864
10.238 13.799
40
5 10
10.873 15.142
10.937 15.244
10.854 15.107
10.789 15.068
4.3. Free vibration analysis of laminated plates 4.3.1. Natural frequency analysis of cross-ply laminates A simply supported cross-ply square laminated composite plate is studied with different modulus ratios (E11/E22). The four-layer plate symmetrically laminated with [0/90/90/0] is used and each layer has the same thickness. The following parameters are used in the computation: L¼10.0 m, t¼1 or 2 m, E22 ¼1.0 GPa, G12 ¼G13 ¼0.6E22, G23 ¼0.5E22, u12 ¼u23 ¼u13 ¼0.25 and r ¼1 kg/m3. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate with 12 12 meshes is used in the analysis. Table 2 shows the non-dimensionalized fundamental frequencies, pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi o ¼ o11 L2 r=ðE22 t2 Þ, of laminates as functions of modulus ratios (E11/E22). Two kinds of length/thickness ratios are considered. Solutions obtained by Reddy [1] and those of Dai et al. [5] are also provided in Table 2 for comparison purpose. It can be seen that the primary frequency increases with the ratios E11/E22 and L/t. The present results are in a very good agreement with those obtained by other researchers. We can also find that the frequencies obtained using present method is a little smaller than DSG. It indicates that the stiffness of model is softened using the edge-based smoothing operation, which compensates nicely the ‘‘overly-stiff’’ behavior of the FEM model as discussed above. The effect of length/thickness ratios on dimensionless fundapffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mental frequencies (o ¼ o11 L2 r=ðE22 t 2 Þ) is shown in Table 3 with a E11/E22 ratio specified of 40. The present results are compared with solutions given by Reddy and Phan [26], as well as Mallikarjuna and Kant [27]. It is observed that good agreements are attained with Reddy and Phan. For thick plate the present results are in between these two solutions of [26,27]. It is worth pointing out that both Reddy and Mallikarjuna use higher-
y ¼ 251
Reddy [1]
Dai et al. [5]
order shear deformation theory; the present method uses the first-order shear deformation theory for this problem.
4.3.2. Free vibration analysis of angle-ply laminated composite plates Two cases of square angle-ply laminated composite plates are further studied in this section. A four-layer angle-ply [y/ y/y/ y] square plate is first considered. The antisymmetric laminate has the same thickness per ply and is simply supported at all edges. The effect of various fiber orientation angles is investigated. The geometry and material properties are given as: L¼10.0 m, t¼1 m, E11 ¼ 40.0 GPa, E22 ¼1.0 GPa, G12 ¼ G13 ¼0.6E22, G23 ¼0.5E22, u12 ¼u23 ¼u13 ¼0.25 and r ¼1 kg/m3. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate with 16 16 meshes is used in the analysis. The exact solutions given by Bert and Chen [28], the finite strip solutions given by Wang and Lin [29], and other finite element solutions given by Mallikarjuna and Kant [27], Reddy [2], are used for comparison. The fundamental frequencies normalized by o ¼ o11 L2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r=ðE22 t2 Þ versus orientation angle y are tabulated in Table 4. It
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
825
Table 3 Effect of aspect ratio L/t on the non-dimensionalized fundamental frequencies o for a simply supported four-layer cross-ply [0/90/90/0] square laminated composite plate. L/t
4
5
10
20
50
100
Present DSG [24] Reddy and Phan [26] Mallikarjuna and Kant [27]
9.406 9.459 9.497 9.258
10.873 10.937 10.989 10.740
15.142 15.244 15.270 15.090
17.605 17.733 17.668 17.637
18.588 18.727 18.606 18.669
18.746 18.885 18.755 18.835
Table 4 Non-dimensional fundamental frequencies o for a simply supported four-layer angle-ply [y/ y/y/ y] square laminated plate with varying fiber orientations.
y
51
151
301
451
Present DSG [24] Bert and Chen [28] Wang and Lin [29] Mallikarjuna and Kant [27] Reddy [2]
14.686 14.743 14.740 14.853 – –
15.540 15.604 15.610 15.759 – –
17.544 17.622 17.630 17.702 17.513 –
18.323 18.405 18.460 18.263 18.320 18.259
Table 5 Non-dimensional fundamental frequencies on for a simply supported ten-layer angle-ply [y/ y/y] square laminated plate with varying aspect L/t. L/t
100 20 10 6.67 5 4 3.33
y ¼ 151
y ¼ 301
domains. The rectangle rule is employed for performing the numerical integration over each smoothing domain. Based on the results from extensive numerical evaluations, the following features of the present method are found. (1) Present ES-FEM uses 3-node triangular element with five degrees of freedom per node. There is no extra sampling point introduced to evaluate the stiffness matrix in present formulation. The formulation is very simple, and its performance is significantly improved by using an edge-based smoothing operation. (2) The solutions obtained using 3-noded triangular elements in the present study are very stable and accurate in comparison to other complex plate elements. (3) Satisfactory accuracy and rate of convergence can be obtained for both thin and thick the laminated composite plates. (4) The present element can be easily extended to quadrangular element and general n-sided (n44) polygonal elements.
y ¼451
Present
Noor and Burton [4]
Present
Noor and Burton [4]
Present
Noor and Burton [4]
0.001323 0.03193 0.1160 0.2303 0.3586 0.4927 0.6290
0.001328 0.03201 0.1162 0.2304 0.3588 0.4934 0.6307
0.001503 0.03611 0.1299 0.2548 0.3925 0.5341 0.6762
0.001510 0.03620 0.1296 0.2532 0.3889 0.5286 0.6692
0.001586 0.03796 0.1353 0.2630 0.4023 0.5445 0.6866
0.001595 0.03808 0.1351 0.2617 0.3993 0.5400 0.6810
can be found that the present results agree well with those compared solutions. The maximum deviation of the present solutions from the exact solutions, which are given by Bert and Chen [28], is only 0.72%. A ten-layer angle-ply [y/ y/]5 laminated composite square plate of length L¼1.0 m is used for studying the effect of length/ thickness ratios. The laminate has the same thickness per layer and is simply supported at all edges. Seven different thickness, t ¼0.01, 0.05, 0.10, 0.15, 0.20, 0.25 and 0.30 m, are considered in this study. The material properties are given as: E11 ¼15.0 GPa, E22 ¼1.0 GPa, G12 ¼G13 ¼ 0.5E22, G23 ¼0.35E22, u12 ¼ u13 ¼0.30, u23 ¼0.49 and r ¼1 kg/m3. The shear correction factors k1 and k2 are taken to be 5/6. The whole plate with 20 20 meshes is used in the analysis. The elasticity solutions obtained by Noor and Burton [4] are used as p the reference solutions. The solutions ffiffiffiffiffiffiffiffiffiffiffiffi ffi normalized by o ¼ o11 t r=E22 are shown in Table 5. It can be found that the present solutions agree well with the reference solutions for all the length/thickness ratios tested here.
5. Conclusion Bending and vibration analysis of laminated composite plates has been studied based on an edge-based smoothing technique for using simple 3-node triangular elements. In the present formulation, the smoothed Galerkin weak form is adopted to build the discretized system equations in edge-based smoothing
Acknowledgments The support of National Science Foundation of China (11002053), National Basic Research Program of China (2010CB328005), Program for Changjiang Scholar and Innovative Research Team in University, and Science Fund of State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body (51075005) are gratefully acknowledged.
References [1] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. CRC press; 2004. [2] Reddy JN. Static free vibration of antisymmetric, angle-ply laminated plates including transverse shear deformation by the finite element method. Journal of Sound and Vibration 1979;66(4):565–76. [3] Lakshminarayana HV, Murthy SS. A shear-flexible triangular finite element model for laminated composite plates. International Journal for Numerical Methods in Engineering 1984;20(4):591–623. [4] Noor AK, Burton WS. Assessment of computational models for multilayered anisotropic plates. Composite Structures 1990;14(3):233–65. [5] Dai KY, Liu GR, Lim KM, Chen XL. A mesh-free method for static and free vibration analysis of shear deformable laminated composite plates. Journal of Sound and Vibration 2004;269(3–5):633–52. [6] Setoodeh AR, Karami G. Static free vibration and buckling analysis of anisotropic thick laminated composite plates on distributed and point elastic supports using a 3-D layer-wise FEM. Engineering Structures 2004;26:211–20. [7] Lee WH, Han SC. Free and forced vibration analysis of laminated composite plates and shells using 9-node assumed strain shell element. Computational Mechanics 2006;39(1):41–58. [8] Liu GR. A generalized gradient smoothing technique and smoothed bilinear form for Galerkin formulation of a wide class of computational methods. International Journal of Computational Methods 2008;5(2):199–236. [9] Chen JS, Wu CT, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin meshfree methods. International Journal for Numerical Methods in Engineering 2001;50:435–66. [10] Liu GR, Zhang GY, Dai KY, Wang YY, Zhong ZH, Li GY, et al. A linearly conforming point interpolation method (LC-PIM) for 2D solid mechanics problems. International Journal of Computational Methods 2005;2:645–65. [11] Liu GR, Li Y, Dai KY, Luan MT, Xue WA. Linearly conforming radial point interpolation method for solid mechanics problems. International Journal of Computational Methods 2006;3(4):401–28. [12] Liu GR, Zhang GY. Upper bound solution to elasticity problems: a unique property of the linearly conforming point interpolation method (LC-PIM). International Journal for Numerical Methods in Engineering 2008;74(7):1128–61.
826
X.Y. Cui et al. / Engineering Analysis with Boundary Elements 35 (2011) 818–826
[13] Liu GR, Zhao X, Dai KY, Zhong ZH, Li GY, Han X. Static and free vibration analysis of laminated composite plates using the conforming radial point interpolation method. Composites Science and Technology 2008;68(2):354–66. [14] Puso MA, Solberg J. A stabilized nodally integrated tetrahedral. International Journal for Numerical Methods in Engineering 2006;67:841–67. [15] Puso MA, Chen JS, Zywicz E, Han X, Elmer W. Meshfree and finite element nodal integration methods. International Journal for Numerical Methods in Engineering 2008;74:416–46. [16] Nagashima T. Node-by-node meshless approach and its applications to structural analyses. International Journal for Numerical Methods in Engineering 1999;46:341–85. [17] Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. Computational Mechanics 2007;39(6):859–77. [18] Liu GR, Nguyen TT, Dai KY, Lam KY. Theoretical aspects of the smoothed finite element method (SFEM). International Journal for Numerical Methods in Engineering 2007;71(8):902–30. [19] Cui XY, Liu GR, Li GY, Zhao X, Nguyen TT, Sun GY. A smoothed finite element method (SFEM) for linear and geometrically nonlinear analysis of plates and shells. CMES: Computer Modeling in Engineering & Science 2008;28(2):109–26. [20] Nguyen-Xuan H, Rabczuk T, Bordas S, Debongnie JF. A smoothed finite element method for plate analysis. Computer Methods in Applied Mechanics and Engineering 2008;197:1184–203. [21] Nguyen-Xuan H, Nguyen-Thoi T. A stabilized smoothed finite element method for free vibration analysis of Mindlin–Reissner plates. Communications in Numerical Methods in Engineering 2009;25:882–906.
[22] Dai KY, Liu GR, Nguyen. TT. An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics. Finite Elements in Analysis and Design 2007;43:847–60. [23] Liu GR, Nguyen TT, Lam. KY. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analysis. Journal of Sound and Vibration 2009;320:1100–30. [24] Bletzinger KU, Bischoff M, Ramm. E. A unified approach for shear-locking-free triangular and rectangular shell finite elements. Computers & Structures 2000;75(3):321–34. [25] Spilker RL, Chou SC, Orringer O. Alternate hybrid-stress elements for analysis of multilayer composite plates. Journal of Composite Materials 1977;11(1):51–70. [26] Reddy JN, Phan ND. Stability and vibration of isotropic, orthotropic, and laminated plates using a higher-order shear deformation theory. Journal of Sound and Vibration 1985;98(2):157–70. [27] Mallikarjuna T Kant. Free vibration of symmetrically laminated plates using a higher order theory with finite element technique. International Journal for Numerical Methods in Engineering 1989;28(8):1875–89. [28] Bert CW, Chen TLC. Effect of shear deformation on vibration of antisymmetric angle-ply laminated rectangular plates. International Journal of Solids and Structures 1978;14(6):465–73. [29] Wang WJ, Lin KC. Free vibration of laminated plates using a finite strip method based on a higher-order plate theory. Computers & Structures 1994;53(6): 1281–9.