- Email: [email protected]
Super convergent shear deformable finite elements for stability analysis of composite beams
Composites: Part B 44 (2013) 100–111
Contents lists available at SciVerse ScienceDirect
Composites: Part B journal homepage: www.elsevier.com/locate/compositesb
Super convergent shear deformable finite elements for stability analysis of composite beams Nam-Il Kim 1, Dong-Ho Choi ⇑ Department of Civil and Environmental Engineering, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea
a r t i c l e
i n f o
Article history: Received 17 November 2011 Received in revised form 26 June 2012 Accepted 6 July 2012 Available online 21 July 2012 Keywords: A. Glass fibers A. Lamina/ply B. Buckling C. Numerical analysis
a b s t r a c t The super convergent finite beam elements are newly presented for the spatially coupled stability analysis of composite beams. For this, the theoretical model applicable to the thin-walled laminated composite I-beams subjected to the axial force is developed. The present element includes the transverse shear and the warping induced shear deformation by using the first-order shear deformation beam theory. The stability equations and force–displacement relationships are derived from the principle of minimum total potential energy. The explicit expressions for the seven displacement parameters are then presented by applying the power series expansions of displacement components to simultaneous ordinary differential equations. Finally, the element stiffness matrix is determined using the force–displacement relationships. In order to demonstrate the accuracy and the superiority of the beam element developed by this study, the numerical solutions are presented and compared with the results obtained from other researchers, the isoparametric beam elements based on the Lagrangian interpolation polynomial, and the detailed three-dimensional analysis results using the shell elements of ABAQUS. The effects of shear deformation, boundary condition, fiber angle change, and modulus ratios on buckling loads are investigated in the analysis. Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved.
1. Introduction The laminated composite I-beams have been used extensively in civil and mechanical engineering as well as in aerospace engineering, where high strength and stiffness, and low weight are of primary importance. The fiber composites have also many desirable characteristics, such as corrosion resistance, magnetic transparency, and excellent fatigue characteristics in the direction of the fibers. However, these laminated composite I-beams might be subjected to an axial force when used in above applications and are very susceptible to flexural–torsional buckling. Therefore, the accurate prediction of their stability limit state is of fundamental importance in the design of composite structures. The theory of the thin-walled beams made of isotropic materials was first conducted by Vlasov [1] and Gjelsvik [2] presented an isotropic beam theory identified with plate segments of the beam. Intensive research works [3–10] have been made over the years to develop the theoretical beam models and obtain analytical solutions for the stability analysis of thin-walled composite beams since the early works of Bauld and Tzeng [3]. Kollár [4] presented ⇑ Corresponding author. Tel.: +82 2 2220 0328. E-mail addresses: [email protected] (N.-I. Kim), [email protected] (D.-H. Choi). 1 Tel.: +82 2 2220 4154.
the stability analysis of axially loaded, thin-walled open section, and orthotropic composite I-beams. In his study, a simplified and approximate solution was derived, in which the effect of shear deformations is included using Föppl’s theorem. Davalos et al. [5] developed the total potential energy equations for the instability of FRP I-beams using nonlinear elastic theory. They formulated the simplified engineering equations for critical flexural–torsional buckling loads. A combined analytical and experimental study of the flexural–torsional buckling of FRP composite cantilever beams with I- and channel-sections was presented by Qiao et al. [6] and Shan and Qiao [7], respectively, based on an energy method. Sherbourne and Kabir [8] presented an analytical view of the lateraltorsional buckling of thin-walled composite and single-layered members, highlighting the transverse shear strain effect. Pandey et al. [9] obtained the beam stiffness coefficients which account for cross-section geometry, the material anisotropy of the section as well as the geometrical characteristics of columns, based on a Vlasov-type linear hypothesis. Cortı´nez and Piovan [10] developed the theoretical model for the stability and dynamic analyses of composite thin-walled beams with open or closed cross-sections and obtained an analytical solution for the simply supported beams. On the other hand, the finite element method has also been widely used because of its versatility and accordingly a large amount of work [11–18] was devoted to the improvement of
1359-8368/$ - see front matter Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compositesb.2012.07.013
101
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
composite finite elements. Vo and Lee [11] presented the flexural– torsional coupled buckling of thin-walled composite beams with arbitrary lay-ups based on the first-order shear deformable beam theory. In their study, the generalized displacements were expressed over each element as a linear combination of the onedimensional Lagrangian interpolation function. Vo and Lee [12] developed the quadratic beam elements with three-node to investigate the coupled buckling and vibration behavior of the shear deformable thin-walled laminated beams. Back and Will [13] developed the first-order shear deformable beam theory and three different types of finite elements, namely, linear, quadratic and cubic elements to solve the governing equations. They used reduced integration to alleviate shear locking. Suresh and Malhotra [14] used the shear deformable plated theory with finite element method to analyze the buckling behavior of the composite box type of configuration. Lin et al. [15] developed a finite element method, in which the element has seven DOFs at each node, to study the stability of thin-walled FRP structural members. Cortı´nez and Piovan [16] presented a theoretical model for the stability analysis of shear deformable composite thin-walled beams with open or closed cross-sections and developed a finite element with twonode and 14 DOFs to solve the governing equations. The generalized beam theory (GBT) formulation to analyze the local and global buckling behavior of FRP composite thin-walled columns with arbitrary open cross-sections was presented by Silva et al. [17] based on the finite element method. Silva and Silvestre [18] developed the GBT formulation incorporating the effects of elastic coupling, warping, cross-section in-plane deformation, and shear deformation and used the finite element formulation to solve the GBT system of differential equilibrium equations. A new theory for the stability analysis of composite beams that included the coupling effects caused by the Poisson effect was presented by Shield and Morey [19] using the Ritz method. They analyzed the I-section and box beams that are middle plane symmetric lay-ups and showed that the inclusion of the anticlastic curvature substantially reduces the predicted buckling loads. Bhaskar and Librescu [20] conducted a study of the flexural buckling of axially compressed thin-walled beams featuring the extention-twist coupling using the energy method based on Trefftz’s criterion (Washizu [21]). They discussed the effects of direct transverse shear and the parasitic bending and transverse shear coupling. Another powerful and effective approach solving the stability problems of composite beam is to develop the stiffness matrix based on solutions of the simultaneous ordinary differential equations. This stiffness matrix approach uses the exact shape functions in the form of converging infinite series. By using this method for any desired set of boundary conditions and assembly of elements, the critical buckling loads can be evaluated. Recently Jun et al. [22] introduced the dynamic stiffness method for analyzing the buckling and free vibration behaviors of axially loaded laminated composite beams having arbitrary lay-up. Abramovich et al. [23] applied an exact element method to calculate the buckling loads and the natural frequencies of nonsymmetric laminated rectangular composite beams. They investigated the influence of longitudinal and transverse restraints, material, and lay-up sequence on the buckling loads and natural frequencies of beam. However, in the above studies, only composite beams with rectangular cross-section were considered in the analysis. The existing literatures reveal that although a large number of studies dealing with the buckling behavior of composite beams have been performed, it should be noted that there was no study reported on the evaluation of stiffness matrix for the stability analysis of the thin-walled laminated composite beams considering the fully coupled shear deformation effects in the literature. It is well known that evaluating exactly the critical buckling loads of the shear deformable thin-walled composite beam having arbitrary
lay-ups is very difficult due to the complexities arising from coupling effects of the flexural and torsional deformation together with shear deformation. Accordingly, the contribution of this study is to present the 14 14 element stiffness matrix of shear deformable thin-walled laminated composite I-beams. The important points are summarized as follows: 1. A general theory is developed for the spatially coupled stability analysis of thin-walled composite beams with arbitrary lay-up considering the transverse shear and the warping induced shear deformation by using the first-order shear deformation beam theory. The Wagner effect which accounts for the change in the effective torsional stiffness due to the axial stresses is considered. 2. Based on the power series expansions of displacement components, the numerical method to exactly evaluate the element stiffness matrix of shear deformable thin-walled composite beam with arbitrary lay-up is presented. 3. For comparison, the finite element procedure based on the Largangian interpolation polynomial is presented. 4. To demonstrate accuracy and validity of this study, the critical buckling loads of I-beams with symmetric and arbitrary layups are evaluated and compared with the results from available references and the finite element solutions using isoparametric beam elements and shell elements from ABAQUS [24]. 2. General theory for shear deformable composite beam 2.1. Variational formulation In this study, two coordinate systems are used as shown in Fig. 1. The first coordinate system is the local coordinate (g, s, x), wherein the g axis is normal to the middle surface of the element and the s axis is tangent to the middle surface and is directed along the contour line of the cross-section. The second coordinate system is the orthogonal Cartesian coordinate (x, y, z), for which y and z axes lie in the plane of the cross-section and the x axis is parallel to the longitudinal axis of the beam; r is the g coordinate and q is the s coordinate of Q; x1 is the rotation of the cross-section about the axis of the beam and w is the angle tangent to the contour makes with y axis. To derive the composite beam theory, the following assumptions are made: 1. The strains are assumed to be small. 2. The beam is linearly elastic and prismatic. 3. The cross-section is assumed to maintain its shape during deformation, so that there is no distortion.
s, v
w Q q
,u
r
P O
Fig. 1. Pictorial definitions of coordinates in thin-walled cross-section.
102
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
4. The transverse shear strains and the warping shear are incorporated. It is assumed that they are uniform over the crosssections. From the study by Kim [25], the displacement fields of any generic point Q on the profile section can be expressed in terms of the seven displacement parameters as follows:
u ¼ U y sin w U z cos w x1 q @u @s w ¼ U x x3 y þ x2 z þ f / gðx3 sin w þ x2 cos w þ fqÞ
v ¼ Uy cos w þ Uz sin w þ x1 r g
dPE ¼
Z
1 2
L
Z A
o
dPG ¼
L
h
PfU 0y dU 0y
þ
U 0z dU 0z
þ ðx
þ
U 0y d
Y
! dt ¼ 0
ð8Þ
ext
t2
Pdt ¼
Z
t2 t1
Z
o
L
h
F 1 dU 0x þ F 2 dðU 0y x3 Þ þ F 3 dðU 0z þ x2 Þ
þ M2 dx02 M 3 dx03 þ M/ df 0 þ Tdðf þ x01 Þ þ Mt dðx01 f Þ þ PfU 0y dU 0y þ U 0z dU 0z þ ðx01 dU 0y þ U 0y dx01 Þzp ðx01 dU 0z þ U 0z dx01 Þyp i þM p x01 dx01 Te dDTe dxdt ¼ 0
ð4Þ
where yp and zp are coordinates of the pole P as shown in Fig. 1 and the following geometric relations are used.
ð9Þ
The constitutive equations [27] of the kth lamina in the local coordinate system of the cross-section are written as
rx rsx
"
k ¼
Q 11
Q 16
Q 16
Q 66
#k
ex csx
ð10Þ
where Q ij are the transformed reduced stiffness coefficients and include the material properties of each lamina. The stress–strain relationship can be simplified by adopting the zero hoop stress assumption (rs = 0) or the zero hoop strain one (es = 0) and the detailed explanation can be found in Ref. [27]. 2.2. Stability equations and force–displacement relationships The stability equations of the present study can be obtained by integrating the derivatives of the varied quantities by parts in Eq. (9), and collecting the coefficients of dUx, dUy, dUz, dx1, dx2, dx3, and df.
E11 U 00x þ E16 U 00y þ E17 U 00z þ ðE15 þ E18 Þx001 þ E12 x002 þ E17 x02 E13 x003 E16 x03 þ E14 f 00 þ ðE18 E15 Þf 0 ¼ 0
ð11aÞ
E16 U 00x þ E66 U 00y þ E67 U 00z þ ðE56 þ E68 Þx001 þ E26 x002 þ E67 x02 E36 x003 E66 x03 þ E46 f 00 þ ðE68 E56 Þf 0 þ PU 00y þ Pzp x001 ¼ 0
0 1 Þzp
x
o ðx01 dU 0z þ U 0z dx01 Þyp g þ Mp x01 dx01 dx
Z
ð3Þ
0 0 1 dU y
G
ð1cÞ
where A is the area of the cross-section. Variation of the potential energy can be expressed by substituting the assumed displacement fields in Eqs. (1a) and (1b) into Eq. (3) as follows:
Z
E
t1
ð2Þ
P 02 ðu þ v 02 ÞdAdx A
þ
Y
d
where L is the length of beam; F1 is the axial force; F2 and F3 are the shear forces in y and z directions, respectively; M2 and M3 are the bending moments about y and z axes, respectively; M/ is the bimoment; T and Mt are the two contributions to the total twisting moment; eox is the axial strain; jy and jz are the biaxial curvatures in y and z directions, respectively; j/ and jsx are the warping curvature and twisting curvature, respectively; coyx and cozx are the transverse shear strains; co/ is the warping shear strain. Now the potential energy due to the initial axial force P acting at the centroid is given as follows:
PG ¼
d t1
Y
ð1bÞ
ðF 1 deox þ F 2 dcoyx þ F 3 dcozx þ M2 djy þ M3 djz þ M/ dj/
þ Tdco/ þ M t djsx Þdx
t2
Substituting Eqs. (2), (4), and (7) into the extended Hamilton’s principle in Eq. (8), the following weak statement is obtained.
L
o
Z
ð1aÞ
where Ux, Uy, and Uz are the rigid body displacements; x1 is the rotation of cross-section; x2 and x3 are the rotations about y and z axes, respectively; f is the warping parameter. The variation of strain energy for the thin-walled composite beams considering the transverse shear and the warping induced shear deformation is expressed as follows [25]:
Z
In order to derive the stability equations and the force–displacement relationships of the laminated composite beam, the extended Hamilton’s principle is used by
ð11bÞ
E17 U 00x þ E67 U 00y þ E77 U 00z þ ðE57 þ E78 Þx001 þ E27 x002 þ E77 x02 E37 x003 E67 x03 þ E47 f 00 þ ðE78 E57 Þf 0 þ PU 00z Pyp x001 ¼ 0
ð11cÞ
ðE15 þ E18 ÞU 00x þ ðE56 þ E68 ÞU 00y þ ðE57 þ E78 ÞU 00z þ ð2E58 þ E55 þ E88 Þx001
y yp ¼ q cos w þ r sin w
ð5aÞ
þ ðE25 þ E28 Þx002 þ ðE57 þ E78 Þx02 ðE35 þ E38 Þx003 ðE56 þ E68 Þx03
z zp ¼ q sin w r cos w
ð5bÞ
þ ðE45 þ E48 Þf 00 ðE55 E88 Þf 0 þ Pzp U 00y Pyp U 00z þ Mp x001 ¼ 0
In Eq. (4), the stress resultant Mp which is known as the Wagner effect is evaluated by
Mp ¼
Z A
P 2 ðq þ r 2 þ g2 þ 2grÞdsdg A
ð6Þ
This property accounts for the change in the effective torsional stiffness due to the bending compressive and tensile stresses, which create a torque as the beam twists during buckling. For the isotropic beams, Eq. (4) is the same as the one in Kim et al. [26]. The variation of the potential energy due to external loads can be written as follows:
dPext ¼
dDTe Te
ð7Þ
where De and Te signify the member end displacement and force vector, respectively.
ð11dÞ
E12 U 00x E17 U 0x þ E26 U 00y E67 U 0y þ E27 U 00z E77 U 0z þ ðE25 þ E28 Þx001 ðE57 þ E78 Þx01 þ E22 x002 E77 x2 E23 x003 þ ðE37 E26 Þx03 þ E67 x3 þ E24 f 00 ðE25 E28 þ E47 Þf 0 þ ðE57 E78 Þf ¼ 0
ð11eÞ
E13 U 00x E16 U 0x þ E36 U 00y E66 U 0y þ E37 U 00z E67 U 0z þ ðE35 þ E38 Þx001 ðE56 þ E68 Þx01 þ E23 x002 þ ðE37 E26 Þx02 E67 x2 E33 x003 þ E66 x3 þ E34 f 00 ðE35 E38 þ E46 Þf 0 þ ðE56 E68 Þf ¼ 0 ð11fÞ E14 U 00x þ ðE15 E18 ÞU 0x þ E46 U 00y þ ðE56 E68 ÞU 0y þ E47 U 00z þ ðE57 E78 ÞU 0z þ ðE45 þ E48 Þx001 þ ðE55 E88 Þx01 þ E24 x002 þ ðE25 E28 þ E47 Þx02 þ ðE57 E78 Þx2 E34 x003 ðE35 E38 þ E46 Þx03 þ ðE68 E56 Þx3 þ E44 f 00 ðE55 2E58 þ E88 Þf ¼ 0
ð11gÞ
103
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
where Eij are the laminate stiffnesses which depend on the crosssection of the composite and given in Ref. [25]. The seven force– displacement relationships are as follows:
The solutions of seven displacement measures are taken as the following infinite power series.
F 1 ¼ E11 U 0x þ E16 U 0y þ E17 U 0z þ ðE15 þ E18 Þx01 þ E17 x2
Ux ¼
þ E12 x02 E16 x3 E13 x03 þ ðE18 E15 Þf þ E14 f 0
ð12aÞ
ð12bÞ
F 3 ¼ E17 U 0x þ E67 U 0y þ ðE77 þ PÞU 0y þ ðE78 þ E57 Pyp Þx01 ð12cÞ
M 1 ¼ ðE15 þ
Pyp ÞU 0z
þ ðE68 þ E56 þ
Pzp ÞU 0y
þ ðE78 þ E57
þ E78 Þx2 þ ðE25 þ E28 Þx02 ðE56 þ E68 Þx3 ðE35 þ E38 Þx þ ðE88 E55 Þf þ ðE45 þ E48 Þf M2 ¼
E12 U 0x
E26 U 0y
þ
þ
E27 U 0z
0
þ ðE25 þ E28 Þx þ E27 x2 ð12eÞ
E23 x02 þ E36 x3 þ E33 x03 þ ðE35 E38 Þf E34 f 0
ð12fÞ
þ E24 x02 E46 x3 E34 x03 þ ðE48 E45 Þf þ E44 f 0
ð12gÞ
The Eq. (11) is the most general form of the coupled buckling for the thin-walled laminated composite beam. If the cross-section is bisymmetric and the material is isotropic, in which the coupling effects are neglected, Eq. (11) can be simplified to the following differential equations [28].
EAU 00x ¼ 0
ð13aÞ
x
þ
¼0
ð13bÞ
GA3 ðU 00z þ x02 Þ þ PU 00z ¼ 0
ð13cÞ
GJ 001 þ GAr ð 001 þ f 0 Þ Mp EI2 002 GA3 ðU 0z þ 2 Þ ¼ 0 EI3 003 þ GA2 ðU 0y 3 Þ ¼ 0
x x x
x
1 X g n xn n¼0
Substituting Eq. (15) into Eq. (11) and shifting the index of power of xn, the stability equations can be compactly expressed in a matrix form as follows:
¼ Zn fan ; anþ1 ; bn ; bnþ1 ; cn ; cnþ1 ; dn ; dnþ1 ; en ; enþ1 ; fn ; fnþ1 ; g n ; g nþ1 gT ð16Þ
ð17Þ
From Eq. (16), the following relation is obtained.
00 1
x ¼0
ð13dÞ
x x
EI/ f 00 GAr ðx01 þ f Þ ¼ 0
¼ Zi Ni fai ; aiþ1 ; bi ; biþ1 ; ci ; ciþ1 ; di ; diþ1 ; ei ; eiþ1 ; fi ; fiþ1 ; g i ; g iþ1 gT ð18Þ
M / ¼ E14 U 0x þ E46 U 0y þ E47 U 0z þ ðE45 þ E48 Þx01 þ E47 x2
PU 00y
n¼0
f ¼
faiþ2 ; biþ2 ; ciþ2 ; diþ2 ; eiþ2 ; fiþ2 ; g iþ2 gT
0 1
M 3 ¼ E13 U 0x E36 U 0y E37 U 0z ðE35 þ E38 Þx01 E37 x2
0 3Þ
n¼0
1 X f n xn ;
a ¼ fao ; a1 ; bo ; b1 ; co ; c1 ; do ; d1 ; eo ; e1 ; fo ; f1 ; g o ; g 1 gT ð12dÞ
þ E22 x02 E26 x3 E23 x03 þ ðE28 E25 Þf þ E24 f 0
GA2 ðU 00y
x3 ¼
where Zn is the 7 14 matrix function. We then define the initial integration constant vector a as
þ ðE55 þ 2E58 þ E88 þ M p Þx01 þ ðE57
0 3
n¼0
1 X x2 ¼ en xn ;
fanþ2 ; bnþ2 ; cnþ2 ; dnþ2 ; enþ2 ; fnþ2 ; g nþ2 gT
þ E77 x2 þ E27 x02 E67 x3 E37 x03 þ ðE78 E57 Þf
E18 ÞU 0x
dn xn ;
n¼0
ð15Þ
þ E67 x2 þ E26 x02 E66 x3 E36 x03 þ ðE68 E56 Þf
þ E47 f 0
n¼0 1 X n¼0
F 2 ¼ E16 U 0x þ ðE66 þ PÞU 0y þ E67 U 0z þ ðE56 þ E68 þ Pzp Þx01 þ E46 f 0
x1 ¼
1 1 1 X X X an xn ; U y ¼ b n xn ; U z ¼ c n xn
The terms for ai+2, bi+2, ci+2, di+2, ei+2, fi+2, and gi+2 converge to zero as i ? 1. The displacement state vector in Eq. (14) is expressed with respect to a by using Eqs. (15) and (18) as follows:
d ¼ Xn a
ð19Þ
where Xn denotes the 14 14 matrix function with the coefficients of Ux, Uy, Uz, x1, x2, x3, and f. In each of these 14 solution sets, the calculation of the coefficients by the recursive relations in Eq. (18) is continued until the contribution of the next coefficient is less than an arbitrarily chosen small number. In this study, above symbolic calculations are performed with the help of the technical computer software Methematica [29]. Next, the initial integration constant vector a can be expressed with respect to 14 nodal displacement components at two ends of the laminated beam. For this, we consider the nodal displacement vector at a and b, which mean the two ends of member, is defined by
ð13eÞ
Ue ¼ hUa ; Ub iT
ð13fÞ
Substituting coordinates of the two ends of member (x = 0, L) into Eq. (19), the nodal displacement vector Ue is expressed as follows:
ð13gÞ
ð20Þ
where EA is the axial rigidity; GA2, and GA3 are the flexural shear rigidities with respect to y and z axes, respectively; GJ and GAr are the torsional rigidity and the warping shear rigidity, respectively; EI2 and EI3 are the flexural rigidities with respect to y and z axes, respectively; EI/ are the warping rigidity.
The displacement state vector consisting of 14 displacement components can be rewritten by using Eqs. (19) and (21) as follows:
3. Exact stiffness matrix of laminated composite beam
d ¼ Xn H1 Ue
3.1. Evaluation of exact displacement function In order to construct the member stiffness matrix of the laminated composite beam, the exact displacement functions are rigorously evaluated based on the seven stability equations derived in previous Chapter. For this, the displacement state vector consisting of 14 displacement parameters is considered as follows:
d¼
D
U x ; U 0x ; U y ; U 0y ; U z ; U 0z ;
0 1;
0 2;
0 0 3; f ; f
x1 ; x x2 ; x x3 ; x
ET
ð14Þ
Ue ¼ Ha
ð21Þ
ð22Þ
It should be noted that XnH1 in Eq. (22) is not an approximate displacement function but an exact one since the displacement state vector satisfies the homogenous form of the coupled equilibrium equations in Eq. (11) exactly. 3.2. Calculation of stiffness matrix In order to calculate the stiffness matrix of the composite member using the displacement function derived in this study, we consider the nodal force vector at two ends a and b defined by
104
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
Fe ¼ hFa ; Fb iT
ð23Þ
The force–deformation relations in Eq. (12) can be rewritten in a matrix form.
f ¼ Sd
ð24Þ
where S is the 7 14 matrix. Substitution of the displacement function in Eq. (22) into Eq. (24) leads to
f ¼ SXn H1 Ue
ð25Þ
Finally, the stiffness matrix of the fully coupled laminated composite beam is evaluated as
Fe ¼ KUe
ð26Þ
where
" K¼
SXn ð0ÞH1 SXn ðLÞH
# ð27Þ
1
The critical buckling loads Pcr of the beam are the values that cause the stiffness matrix for the element to become singular as Eq. (27). A search procedure is employed to find these values up to the desired accuracy. In this study, the Regular-Falsi method [30] is applied to ensure that none of the buckling loads is missed.
det jKj ¼ 0
ð28Þ
It is noteworthy that the stiffness matrix in Eq. (27) is formed by the shape functions which are exact solutions of the governing equations. Therefore, the thin-walled composite beam based on the stiffness matrix developed by this study eliminates discretization errors and is free from the shear and membrane locking. 4. Finite element formulation
The element displacement vector Uie and the force vector Fie for the isoparametric beam element are defined as
Uie ¼ ½U 1 ; U 2 n
U ¼
½U nx ; U ny ; U nz ;
ð30aÞ n 1;
n 2;
x x x
n n T 3; f ;
n ¼ 1; 2
ð30cÞ ð30dÞ
where Un and Fn are the nodal displacement and force vectors of the isoparametric beam element, respectively. Substituting the shape functions and the cross-sectional properties into Eq. (9) and integrating along the element length, the total potential energy of the finite composite beam element is obtained in a matrix form as
Y
¼
1 T U ðKie þ Kig ÞUie UTie Fie 2 ie
ð31Þ
where Kie and Kig are the element elastic and geometric stiffness matrices in local coordinate, respectively. When equal interpolation of all seven displacement parameters and full integration is used to evaluate the stiffness matrices Kie and Kig, the beam element is known to behavior very stiff. This behavior is due to the instability of the element to represent a constant state of the transverse shear strains coyx ; cozx and the warping shear strain co/ , and the phenomenon is known as the shear locking. To avoid locking, the shear strains are approximated, while using equal interpolation of the seven displacement parameters as a constant by using the reduced integration of the shear stiffnesses. Therefore, in the present work, the reduced integration scheme in which the quadrature rule is taken to be one order lower than the normal one is employed. For the entire structure, the assembly of element stiffness matrix leads to the equilibrium matrix equation in a global coordinate system based on the coordinate transformation.
For comparison, the finite element formulation based on the shear deformable isoparametric beam element with seven DOFs per a node is presented. The present finite beam element uses the same shape functions for all translational, rotational, and warping displacements. The displacements can be expressed as follows:
l ð1 þ r Þ 2 2 X Ui ¼ Rn ðr ÞU ni
x¼
ð30bÞ
Fie ¼ ½F 1 ; F 2 h iT Fn ¼ F n1 ; F n2 ; F n3 ; M n1 ; Mn2 ; M n3 ; Mn/
1
ð29aÞ i ¼ x; y; z
ð29bÞ
j ¼ 1; 2; 3
ð29cÞ
n¼1
xj ¼
2 X
Rn ðr Þxnj
n¼1
f ¼
2 X Rn ðr Þf n
ð29dÞ
n¼1
where U ni ; xni , and fn are the translational and rotational displacements and the warping parameter at node n, respectively; Rn is the Lagrangian interpolation function whose detailed expression is presented in Ref. [31]; r⁄ is a natural coordinate that varies from 1 to +1.
2
Fig. 2. Shape of cross-section under consideration.
Table 1 Critical buckling loads of the SS bisymmetric I-beams (106 N), (L = 6 m). Stacking sequence
[0]4 [45/ 45]S [0/90]S
Machado and Cortı´nez [32]
This study Without shear
With shear
rs = 0
rs = 0
42.69 4.45 22.90
33.18 4.44 19.85
Vo and Lee [11]
Without shear
With shear
es = 0
rs = 0
rs = 0
rs = 0
33.34 12.99 19.87
42.11 4.45 22.57
33.18 4.44 19.84
30.38 4.41 20.63
105
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
ðKE þ kKG ÞU ¼ F ¼ 0
ð32Þ
Table 2 Critical buckling loads of the CF bisymmetric I-beams (N), (L/h = 20).
where KE and KG are the global elastic and geometric stiffness matrices, respectively; F is the global nodal force vector; k is the load parameter under the assumption of proportional loading.
Stacking sequence
[0]16 [15/15]4S [30/30]4S [45/45]4S [60/60]4S [75/75]4S [0/90]4S
5. Numerical examples In order to illustrate accuracy and validity of the current study, numerical analyses have been performed for the coupled stability analysis of composite I-beams with the bi-symmetric and monosymmetric cross-sections. The terms ‘‘bisymmetric’’ and ‘‘monosymmetric’’ are defined in the geometry of the cross-section not the material point of view. The solutions obtained from this study are compared with those from other researchers and the finite element solutions using the isoparametric beam elements and the shell elements by ABAQUS.
Without shear rs = 0
With shear
rs = 0
es = 0
5755.2 5199.8 3861.0 2672.7 2114.7 1948.3 3857.8
5737.5 5187.7 3856.0 2670.6 2113.2 1946.6 3849.8
5859.1 5413.1 4325.2 3160.2 2370.6 2031.3 3910.8
Vo and Lee [11]
ABAQUS
5741.5 5189.0 3854.5 2668.4 2111.3 1945.1 3829.8
5719.2 5175.3 3853.3 2671.4 2113.4 1945.8 3843.4
5.1. Bisymmetric I-beams with symmetric lay-up The example considered is the simply supported (SS) bisymmetric I-beam, as shown in Fig. 2, whose geometric properties
0.5
7
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
0.4
6
Dimensionless buckling load, Pcr*
Dimensionless buckling load, Pcr*
This study
0.3
0.2
0.1
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
5
4
3
2
1
0.0
0 0
15
30
45
60
75
90
0
30
45
60
Fiber angle, (degree)
(a) CF beams
(c) CS beams
75
90
7
7
6
5
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
6
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
Dimensionless buckling load, Pcr*
Dimensionless buckling load, Pcr*
15
Fiber angle, (degree)
4
3
2
5
4
3
2
1
1
0
0 0
15
30
45
60
75
90
0
15
30
45
60
Fiber angle, (degree)
Fiber angle, (degree)
(b) SS beams
(d) CC beams
Fig. 3. Critical buckling loads of the bisymmetric I-beams with respect to the fiber angle change.
75
90
106
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
are as follows: L = 6000 mm, b1 = b2 = 300 mm, h = 600 mm, and t = 30 mm. The cross-section is made with graphite-epoxy (AS4/ 3501) material and its material constants are as follows: E1 = 144 GPa, E2 = E3 = 9.65 GPa, G12 = G13 = 4.14 GPa, G23 = 3.45 GPa, m12 = m13 = 0.3, m23 = 0.5. In which subscripts ‘1’ and ‘2’ correspond to directions parallel and perpendicular to fibers, respectively. All constituent flanges and web are made of four plies and are assumed to be symmetrically laminated with respect to its middle plane. The critical buckling loads obtained from the present analysis using a single element are compared with the finite element results of Vo and Lee [11] and Machado and Cortı´nez [32] in Table 1. It is observed that the present solutions applying the zero hoop stress assumption (rs = 0) are in good agreement with results in Refs. [11,32]. Next, the bisymmetric I-beams made with glass–epoxy material with the following elastic constants: E1 = 53.78 GPa, E2 = E3 =
17.93 GPa, G12 = G13 = 8.96 GPa, G23 = 3.45 GPa, m12 = m13 = 0.25, m23 = 0.34 is considered. The total 16 layers with equal thickness are considered in two flanges and web and the sectional properties are as follows: b1 = b2 = 25 mm, h = 50 mm, and t = 2.08 mm. For clamped-free (CF) beams with L/h = 20, the buckling loads from this study with and without shear deformation effect are presented in Table 2. For comparison, the finite element solutions by Vo and Lee [11] and ABAQUS are given in Table 2. For ABAQUS calculation, a total of 240 nine node shell elements (S9R5) (six through the cross-section) are used to obtain the results. It can be found from Table 2 that the correlation of the present results considering the shear based on rs = 0 assumption and the solutions from ABAQUS and Vo and Lee [11] is seen to be excellent for all lay-ups considered. It can also be found that the analysis based on es = 0 assumption is seen to overestimate the buckling loads by a maximum 18.3% at w = 45°.
100
1.6
90
70
Dimensionless buckling load, Pcr*
Effect of shear, (%)
1.4
CF SS CS CC
80
60 50 40 30 20
E1/E2 1 3 5 10
1.2 1.0 0.8 0.6 0.4 0.2
10
0.0
0 0
15
30
45
60
75
90
0
15
30
45
60
Fiber angle, (degree)
Fiber angle, (degree)
(a) L h = 5
(a) Critical buckling load
75
90
75
90
0.20
6
0.19
Effect of shear, (%)
4
Dimensionless buckling load, Pcr*
CF SS CS CC
5
3
2
1
E1/E2 1 3 5 10
0.18 0.17 0.16 0.15 0.14 0.13 0.12
0 0
15
30
45
60
75
90
0
15
30
45
60
Fiber angle, (degree)
Fiber angle, (degree)
(b) L h = 20
(b) Region of interest
Fig. 4. Effect of shear for the bisymmetric I-beams with respect to the fiber angle change.
Fig. 5. Critical buckling loads for the flexural mode in y axis of the CF bisymmetric Ibeams with various modulus ratios, (L/h = 15).
107
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
To investigate the effects of boundary condition and shear deformation on the critical buckling loads of the laminated composite beam, the dimensionless buckling loads for four boundary conditions are plotted in Fig. 3 with respect to the fiber angle change. The dimensionless buckling load parameter 2 P cr ¼ P cr L2 =ðAE2 h Þ is introduced for simplicity and the boundary conditions under consideration are: the clamped–free (CF), simply–simply (SS), clamped–simply (CS), and clamped–clamped (CC) boundary conditions. It is seen that the buckling loads of beams without shear are the highest at w = 0°, both of short and long beams for all boundary conditions. Whereas for CC beam with shear, it is the highest around w = 14° for L/h = 5 beam with shear. Fig. 4 shows the variation of the effect of shear on the buckling loads for L/h = 5 and 20 beams. It is observed that the shear effects are the largest at w = 0° and the smallest at w = 50° for both of beams regardless of the boundary conditions. Moreover the maximum shear effect, as can be seen in Fig. 4a, is 5% for CC beam with L/h = 20 which may be relatively slender. Therefore the shear effect 5.0
Dimensionless buckling load, Pcr*
4.5 E1/E2
4.0
1 3 5 10
3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0
15
30
45
60
75
90
Fiber angle, (degree)
(a) Critical buckling load
should be included in the flexural–torsional stability analysis of the laminated composite beam. The effects of the modulus ratio E1/E2 on the critical buckling loads of the CF beams with L/h = 15 are presented in Figs. 5–7. It is seen that the maximum buckling loads of the isotropic beam (E1/E2 = 1) are obtained at w = 45° for the flexural modes in y and z axes as shown in Figs. 5 and 6. While those are exhibited at w = 0° for the anisotropic beams. It is also seen that the critical buckling load increases as the modulus ratio increases in fiber angles less than 60°. Whereas it is reverse in region beyond w = 60°. For the case of torsional mode as shown in Fig. 7, the torsional buckling loads are maximum and minimum at w = 90° for the isotropic and anisotropic beams, respectively. Fig. 8 shows the effect of shear on the flexural and torsional buckling loads. It is evident that the shear effect on the flexural buckling mode about the strong axis is larger than that about the weak axis. It is interesting to observe that for the isotropic beam, the maximum shear effect takes place at w = 45° for the two flexural modes and torsional mode. On the other hand, for anisotropic beams, it occurs at w = 0° for all buckling modes. It is also observed that the shear effect increases as E1/E2 increases in fiber angles w 6 25° and it decreases with increase of E1/E2 in w > 25°. Finally, the convergence study for the isoparametric beam element using the Lagrangian interpolation polynomials and the beam element developed by the present stiffness matrix method is performed and given in Table 3. In case of the stiffness matrix method, the system matrix can be obtained by the usual assembling process from the conditions of equilibrium and compatibility at the nodes. It can be found from Table 3 that the solutions developed by this study using one element are the same as results using three elements, and are in an excellent agreement with those using at least 60 isoparametric beam elements. It should be noted that both the present stiffness matrix method and the isoparametric beam formulation use the same total potential energy. The difference of two methods is that the current element is based on the shape functions which satisfy the homogeneous forms of the equilibrium equations exactly. Therefore it is possible to obtain the solutions though only a minimum number of beam elements are used. However, the solutions obtained from the isoparametric beam elements using the Lagrangian interpolation function which satisfies only displacement continuity at nodal point are approximate. Resultantly, as a number of isoparametric beam element
0.75 2.0 1.8
E1/E2 1 3 5 10
0.65
Dimensionless buckling load, Pcr*
Dimensionless buckling load, Pcr*
0.70
0.60
0.55
0.50
0.45
E1/E2 1 3 5 10
1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2
0.40 0
15
30
45
60
75
90
Fiber angle, (degree)
(b) Region of interest Fig. 6. Critical buckling loads for the flexural mode in z axis of the CF bisymmetric Ibeams with various modulus ratios, (L/h = 15).
0.0 0
15
30
45
60
75
90
Fiber angle, (degree) Fig. 7. Critical buckling loads for the torsional mode of the CF bisymmetric I-beams with various modulus ratios, (L/h = 15).
108
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
2.0
14
1.8 12 1 3 5 10
1.4 1.2
1 3 5 10
10
Effect of shear, (%)
Effect of shear, (%)
E1/E2
E1/E2
1.6
1.0 0.8 0.6
8
6
4
0.4 2 0.2 0
0.0 0
15
30
45
60
75
90
0
15
30
45
60
75
Fiber angle, (degree)
Fiber angle, (degree)
(a) Flexural mode in y axis
(b) Flexural mode in z axis
90
2.0 1.8 E1/E2
Effect of shear, (%)
1.6
1 3 5 10
1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
15
30
45
60
75
90
Fiber angle, (degree)
(c) Torsional mode Fig. 8. Effect of shear of the CF bisymmetric I-beams with various modulus ratios, (L/h = 15).
Table 3 Convergence of the finite beam element and the present solutions for the buckling loads of the CF bisymmetric I-beam with [0/45/90/45]2S lay-up (N), (L/h = 20). Mode
1 2 3
No. of isoparametric beam elements
No. of elements by this study
ABAQUS
5
10
30
60
1
2
3
3379.9 11732.6 11818.3
3338.2 11589.1 11776.7
3326.0 11547.1 11764.6
3324.9 11543.2 11763.5
3324.5 11541.9 11763.1
3324.5 11541.9 11763.1
3324.5 11541.9 11763.1
used in the beam structure increases, its solutions converge to those from the present study using a minimum number of beam elements.
5.2. Monosymmetric I-beams with symmetric and arbitrary lay-ups In this example, the monosymmetric I-beams with symmetric and arbitrary lay-ups are considered. The material is the glass–
3321.3 11350.0 11478.0
epoxy which is used in Section 5.1 and the geometric properties are b1 = 15 mm, b2 = 25 mm, h = 50 mm, and t = 2.08 mm. First, for the CF monosymmetric I-beam with symmetric lay-up, the coupled buckling loads are compared with those from ABAQUS analysis in Table 4 for various stacking sequences. In Table 4, the term ‘‘NW’’ means the buckling load from neglecting warping. It can be found from Table 4 that the good correlation of the current study including rs = 0 assumption and considering warping effect with the ABAQUS analysis is obtained over the whole range of fiber
109
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111 Table 4 Critical buckling loads of the CF monosymmetric I-beams (N), (L/h = 20). Stacking sequence
This study
[0]16 [15/15]4S [30/30]4S [45/45]4S [60/60]4S [75/75]4S [0/90]4S
1.2 ABAQUS
Without shear rs = 0
With shear
rs = 0
rs = 0(NW)
es = 0
3005.9 2854.6 2201.0 1562.4 1241.7 1134.6 2114.0
2992.5 2809.4 2199.0 1561.5 1240.9 1133.8 2111.1
3492.7 3157.6 2346.5 1625.0 1285.8 1184.4 2343.2
3046.4 2917.2 2447.5 1834.9 1386.5 1181.2 2141.0
1.0 0.8
2812.7 2664.9 2129.2 1534.8 1226.6 1122.4 2055.7
0.6
U U U
ω ω ω f
0.4 0.2 0.0
Table 5 Critical buckling loads of the CF monosymmetric I-beam with [0/30/60/90] lay-up (N), (L/h = 15). Mode
1 2 3
This study
ABAQUS
Without shear
With shear
rs = 0
rs = 0
es = 0
3171.1 10951.8 17321.8
3135.9 10900.1 16851.8
3271.9 11176.9 17688.9
-0.2 -0.4 0.00
0.75
1.00
0.75
1.00
(a) Mode 1 2921.1 11189.0 16391.0
Table 6 Critical buckling loads of the CC monosymmetric I-beam with [0/30/60/90] lay-up (N), (L/h = 15).
1 2 3
0.50
x/ L
1.2 1.0
Mode
0.25
0.8 0.6
Without shear
With shear
rs = 0
rs = 0
es = 0
29792.7 57543.3 93819.4
28842.0 54631.8 89936.0
30017.1 56952.8 94770.5
U U U
ω ω ω f
0.4 0.2 0.0
angle. The maximum deviation with es = 0 assumption is about 17.5% at 45° of fiber angle and the effect of warping is the largest at w = 0° with 16.7% and the smallest at w = 60° with 3.6%. It can also be noticed that the present approach exhibits better results than that of classical beam assumption which does not consider shear effect. The purpose of our final example is to evaluate exactly the coupled buckling load of the monosymmetric I-beam with arbitrary lay-up. We consider the monosymmetric I-beam with the same dimensions used in the previous example, but four layers of [0/ 30/60/90] with equal ply thicknesses are used for both flanges and web. Since the lay-up is nonsymmetric with respect to the middle plane of the plate wall, the extension-shear (A16 and A26), bending–twisting (D16 and D26), and extension–twisting and bending–shearing (B16 and B26) couplings have non-zero values. Resultantly all coupling laminate stiffnesses in Eqs. (11) and (12) are non-zero except for E12,E67, and E78. For CF and CC beams with L/ h = 15, the lowest three coupled buckling loads are presented in Tables 5 and 6, respectively. It is seen that the results from the present rs = 0 assumption are in close with ABAQUS solutions for CF beam. The shear effect reduces the coupled buckling loa ds about 2.8% in the 3rd mode for CF beam and 5.3% in the 2nd mode for CC beam. It should be noted that the present solutions obtained from one element are exact for the higher buckling modes as well as the lower ones, while the considerable numbers of finite beam elements based on approximate shape functions may be required to achieve sufficiently accurate results in the higher modes. Furthermore for beams with L/h = 10, the buckling mode shapes corresponding to the lowest two buckling loads are illustrated in
-0.2 -0.4 0.00
0.25
0.50
x/ L
(b) Mode 2 Fig. 9. Coupled buckling mode shapes of CF monosymmetric I-beam with [0/30/60/ 90] lay-up, (L/h = 10).
Figs. 9 and 10 for CF and CC beams, respectively. From Figs. 9 and 10, the flexural–torsional coupled buckling modes and the flexural–torsional-shearing coupled ones can be observed for CF and CC beams, respectively.
6. Conclusion The super convergent finite composite beam elements considering shear effects has been developed for the coupled stability analysis. The stability equations and force–displacement relationships are derived by using the extended Hamilton’s principle and the element stiffness matrix is rigorously evaluated based on the power series expansions of displacement components. For comparison, the finite element model using the Lagrangian interpolation polynomial is presented. An advantage, which is often overlooked but may be more important, is that the present method can provide
110
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
1.2 U U U
1.0
ω ω ω
0.8 0.6
f
0.4 0.2 0.0 -0.2 -0.4 -0.6 0.00
0.25
0.50
0.75
1.00
x/ L
(a) Mode 1
(2) For both of the shorter and longer span beams with bisymmetric cross-section, the minimum shear effect is reached at w = 55° regardless of the boundary conditions. (3) The torsional buckling loads of beams with high modulus ratio are larger than those with small modulus ration whole fiber angles except w = 90°. But it is not always a case of the flexural buckling behavior. (4) For flexural and torsional bucking behavior, the shear effect with high modulus ratio is not always larger than that with small modulus ratio. (5) The present numerical procedure provides a refined method for not only the evaluation of the stiffness matrix of shear deformable laminated composite beam but also general solutions of simultaneous ordinary differential equations of the higher order. This thin-walled composite beam element eliminates discretization errors and is free from locking phenomenon since the displacement state vector satisfies the homogenous form of the simultaneous ordinary differential equations. Furthermore this study is capable of predicting an infinite number of buckling loads of composite beams by means of a finite number of coordinates.
1.2 1.0 0.8 0.6 0.4
U U U
Acknowledgements This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2012R1A1A2007054).
ω ω ω f
0.2
References 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 -1.2 0.00
0.25
0.50
0.75
1.00
x/ L
(b) Mode 2 Fig. 10. Coupled buckling mode shapes of CC monosymmetric I-beam with [0/30/ 60/90] lay-up, (L/h = 10).
benchmark results when compared with other results obtained by the finite element or other approximate methods. Through the numerical examples, the present study is validated comparing the coupled buckling responses with results from other researchers and the finite element solutions using the isoparametric beam elements and the shell elements of ABAQUS. Good correlation is achieved for various beams with different cross-sections, lamination schemes, boundary conditions considered in this paper. The conclusions drawn from this study are as follows: (1) The flexural buckling behavior is affected by the shear and the span-to-height ratio as the fiber angle varies: for shorter span beam with CC boundary condition and with shear effect, the maximum buckling load is obtained at other fiber angle except w = 0°. Whereas for longer span beams with and without shear, and for shorter beam without shear, the maximum buckling loads occur at w = 0°.
[1] Vlasov VZ. Thin walled elastic beams. 2nd ed. Jerusalem: Israel Program for Scientific Transactions; 1961. [2] Gjelsvik A. The theory of thin-walled bars. Wiley; 1981. [3] Bauld NR, Tzeng L. A Vlasov theory for fiber-reinforced beams with thin-walled open cross sections. Int J Solids Struct 1984;20:277–97. [4] Kollá LP. Flexural–torsional buckling of open section composite columns with shear deformation. Int J Solids Struct 2001;38:7525–41. [5] Davalos JF, Qiao P, Salim HA. Flexural–torsional buckling of pultruded fiber reinforced plastic composite I-beams: experimental and analytical evaluations. Compos Strcut 1997;38:241–50. [6] Qiao P, Zou G, Davalos JF. Flexural–torsional buckling of fiber-reinforced plastic composite cantilever I-beams. Compos Strcut 2003;60:205–17. [7] Shan L, Qiao P. Flexural–torsional buckling of fiber-reinforced plastic composite open channel beams. Compos Strcut 2005;68:211–24. [8] Sherbourne AN, Kabir MZ. Shear strain effects in lateral stability of thin-walled fibrous composite beams. J Eng Mech 1995;121:640–7. [9] Pandey MD, Kabir MZ, Sherbourne AN. Flexural–torsional stability of thinwalled composite I-section beams. Compos Eng 1995;5:321–42. [10] Cortı´nez V, Piovan MT. Vibration and buckling of composite thin-walled beams with shear deformability. J Sound Vib 2002;258:701–23. [11] Vo TP, Lee J. On sixfold coupled buckling of thin-walled composite beams. Compos Struct 2009;90:295–303. [12] Vo TP, Lee J. Flexural–torsional coupled vibration and buckling of thin-walled open section composite beams using shear-deformable beam theory. Int J Mech Sci 2009;51:631–41. [13] Back SY, Will KM. Shear-flexible thin-walled element for composite I-beams. Eng Struct 2008;30:1447–58. [14] Suresh R, Malhotra SK. Some studies on buckling of laminated composite thin walled box beams. Compos Struct 1998;40:267–75. [15] Lin ZM, Polyzois D, Shah A. Stability of thin-walled pultruded structural members by finite element method. Thin-Walled Struct 1996;24:1–18. [16] Cortı´nez V, Piovan MT. Stability of composite thin-walled beams with shear deformability. Comput Struct 2006;84:978–90. [17] Silva NMF, Silvestre N, Camotim D. GBT formulation to analyze the buckling behavior of FRP composite open-section thin-walled columns. Compos Strcut 2010;93:79–92. [18] Silva NMF, Silvestre N. On the influence of material couplings on the linear and buckling behavior of I-section composite columns. Int J Strcut Stab Dyn 2007;7:243–72. [19] Shield CK, Morey TA. Kinematic theory for buckling of open and closed section thin-walled composite beams. J Eng Mech 1997;123:1070–81. [20] Bhaskar K, Librescu L. Buckling under axial compression of thin-walled composite beams exhibiting extension-twist coupling. Compos Struct 1995;31:203–12.
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111 [21] Washizu K. Variational method in elasticity and plasticity. Pergammon Press; 1975. [22] Jun L, Hongxing H, Rongying S. Dynamic stiffness analysis for free vibrations of axially loaded laminated composite beams. Compos Struct 2008;84: 87–98. [23] Abramovich H, Eisenberger M, Shulepov O. Vibration and buckling of cross-ply nonsymmetric laminated composite beams. AIAA J 1996;34:1064–9. [24] ABAQUS. Standard user’s manual, Ver. 6.1. Hibbit, Kalsson & Sorensen Inc; 2003. [25] Kim NI. Shear deformable doubly- and mono-symmetric composite I-beams. Int J Mech Sci 2011;53:31–41.
111
[26] Kim MY, Chang SP, Kim SB. Spatial stability and free vibration of shear flexible thin-walled elastic beams. I: analytical approach. Int J Numer Meth Eng 1994;37:4097–115. [27] Jones RM. Mechanics of composite material. Taylor and Francis; 1999. [28] Kim MY, Kim NI, Yun HT. Exact dynamic and static stiffness matrices of shear deformable thin-walled beam-columns. J Sound Vib 2003;267:29–55. [29] Wolfram Mathematica 7. Wolfram Research Inc., IL; 2009. [30] Wendroff B. Theoretical numerical analysis. NY: Academic Press; 1966. [31] Bathe KJ. Finite element procedures. Englewood Cliffs, NJ: Prentice-Hall; 1996. [32] Machado SP, Cortı´nez VH. Non-linear model for stability of thin-walled composite beams with shear deformation. Thin-Walled Struct 2005;43:1615–45.
Contents lists available at SciVerse ScienceDirect
Composites: Part B journal homepage: www.elsevier.com/locate/compositesb
Super convergent shear deformable finite elements for stability analysis of composite beams Nam-Il Kim 1, Dong-Ho Choi ⇑ Department of Civil and Environmental Engineering, Hanyang University, 17 Haengdang-dong, Seongdong-gu, Seoul 133-791, South Korea
a r t i c l e
i n f o
Article history: Received 17 November 2011 Received in revised form 26 June 2012 Accepted 6 July 2012 Available online 21 July 2012 Keywords: A. Glass fibers A. Lamina/ply B. Buckling C. Numerical analysis
a b s t r a c t The super convergent finite beam elements are newly presented for the spatially coupled stability analysis of composite beams. For this, the theoretical model applicable to the thin-walled laminated composite I-beams subjected to the axial force is developed. The present element includes the transverse shear and the warping induced shear deformation by using the first-order shear deformation beam theory. The stability equations and force–displacement relationships are derived from the principle of minimum total potential energy. The explicit expressions for the seven displacement parameters are then presented by applying the power series expansions of displacement components to simultaneous ordinary differential equations. Finally, the element stiffness matrix is determined using the force–displacement relationships. In order to demonstrate the accuracy and the superiority of the beam element developed by this study, the numerical solutions are presented and compared with the results obtained from other researchers, the isoparametric beam elements based on the Lagrangian interpolation polynomial, and the detailed three-dimensional analysis results using the shell elements of ABAQUS. The effects of shear deformation, boundary condition, fiber angle change, and modulus ratios on buckling loads are investigated in the analysis. Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved.
1. Introduction The laminated composite I-beams have been used extensively in civil and mechanical engineering as well as in aerospace engineering, where high strength and stiffness, and low weight are of primary importance. The fiber composites have also many desirable characteristics, such as corrosion resistance, magnetic transparency, and excellent fatigue characteristics in the direction of the fibers. However, these laminated composite I-beams might be subjected to an axial force when used in above applications and are very susceptible to flexural–torsional buckling. Therefore, the accurate prediction of their stability limit state is of fundamental importance in the design of composite structures. The theory of the thin-walled beams made of isotropic materials was first conducted by Vlasov [1] and Gjelsvik [2] presented an isotropic beam theory identified with plate segments of the beam. Intensive research works [3–10] have been made over the years to develop the theoretical beam models and obtain analytical solutions for the stability analysis of thin-walled composite beams since the early works of Bauld and Tzeng [3]. Kollár [4] presented ⇑ Corresponding author. Tel.: +82 2 2220 0328. E-mail addresses: [email protected] (N.-I. Kim), [email protected] (D.-H. Choi). 1 Tel.: +82 2 2220 4154.
the stability analysis of axially loaded, thin-walled open section, and orthotropic composite I-beams. In his study, a simplified and approximate solution was derived, in which the effect of shear deformations is included using Föppl’s theorem. Davalos et al. [5] developed the total potential energy equations for the instability of FRP I-beams using nonlinear elastic theory. They formulated the simplified engineering equations for critical flexural–torsional buckling loads. A combined analytical and experimental study of the flexural–torsional buckling of FRP composite cantilever beams with I- and channel-sections was presented by Qiao et al. [6] and Shan and Qiao [7], respectively, based on an energy method. Sherbourne and Kabir [8] presented an analytical view of the lateraltorsional buckling of thin-walled composite and single-layered members, highlighting the transverse shear strain effect. Pandey et al. [9] obtained the beam stiffness coefficients which account for cross-section geometry, the material anisotropy of the section as well as the geometrical characteristics of columns, based on a Vlasov-type linear hypothesis. Cortı´nez and Piovan [10] developed the theoretical model for the stability and dynamic analyses of composite thin-walled beams with open or closed cross-sections and obtained an analytical solution for the simply supported beams. On the other hand, the finite element method has also been widely used because of its versatility and accordingly a large amount of work [11–18] was devoted to the improvement of
1359-8368/$ - see front matter Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compositesb.2012.07.013
101
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
composite finite elements. Vo and Lee [11] presented the flexural– torsional coupled buckling of thin-walled composite beams with arbitrary lay-ups based on the first-order shear deformable beam theory. In their study, the generalized displacements were expressed over each element as a linear combination of the onedimensional Lagrangian interpolation function. Vo and Lee [12] developed the quadratic beam elements with three-node to investigate the coupled buckling and vibration behavior of the shear deformable thin-walled laminated beams. Back and Will [13] developed the first-order shear deformable beam theory and three different types of finite elements, namely, linear, quadratic and cubic elements to solve the governing equations. They used reduced integration to alleviate shear locking. Suresh and Malhotra [14] used the shear deformable plated theory with finite element method to analyze the buckling behavior of the composite box type of configuration. Lin et al. [15] developed a finite element method, in which the element has seven DOFs at each node, to study the stability of thin-walled FRP structural members. Cortı´nez and Piovan [16] presented a theoretical model for the stability analysis of shear deformable composite thin-walled beams with open or closed cross-sections and developed a finite element with twonode and 14 DOFs to solve the governing equations. The generalized beam theory (GBT) formulation to analyze the local and global buckling behavior of FRP composite thin-walled columns with arbitrary open cross-sections was presented by Silva et al. [17] based on the finite element method. Silva and Silvestre [18] developed the GBT formulation incorporating the effects of elastic coupling, warping, cross-section in-plane deformation, and shear deformation and used the finite element formulation to solve the GBT system of differential equilibrium equations. A new theory for the stability analysis of composite beams that included the coupling effects caused by the Poisson effect was presented by Shield and Morey [19] using the Ritz method. They analyzed the I-section and box beams that are middle plane symmetric lay-ups and showed that the inclusion of the anticlastic curvature substantially reduces the predicted buckling loads. Bhaskar and Librescu [20] conducted a study of the flexural buckling of axially compressed thin-walled beams featuring the extention-twist coupling using the energy method based on Trefftz’s criterion (Washizu [21]). They discussed the effects of direct transverse shear and the parasitic bending and transverse shear coupling. Another powerful and effective approach solving the stability problems of composite beam is to develop the stiffness matrix based on solutions of the simultaneous ordinary differential equations. This stiffness matrix approach uses the exact shape functions in the form of converging infinite series. By using this method for any desired set of boundary conditions and assembly of elements, the critical buckling loads can be evaluated. Recently Jun et al. [22] introduced the dynamic stiffness method for analyzing the buckling and free vibration behaviors of axially loaded laminated composite beams having arbitrary lay-up. Abramovich et al. [23] applied an exact element method to calculate the buckling loads and the natural frequencies of nonsymmetric laminated rectangular composite beams. They investigated the influence of longitudinal and transverse restraints, material, and lay-up sequence on the buckling loads and natural frequencies of beam. However, in the above studies, only composite beams with rectangular cross-section were considered in the analysis. The existing literatures reveal that although a large number of studies dealing with the buckling behavior of composite beams have been performed, it should be noted that there was no study reported on the evaluation of stiffness matrix for the stability analysis of the thin-walled laminated composite beams considering the fully coupled shear deformation effects in the literature. It is well known that evaluating exactly the critical buckling loads of the shear deformable thin-walled composite beam having arbitrary
lay-ups is very difficult due to the complexities arising from coupling effects of the flexural and torsional deformation together with shear deformation. Accordingly, the contribution of this study is to present the 14 14 element stiffness matrix of shear deformable thin-walled laminated composite I-beams. The important points are summarized as follows: 1. A general theory is developed for the spatially coupled stability analysis of thin-walled composite beams with arbitrary lay-up considering the transverse shear and the warping induced shear deformation by using the first-order shear deformation beam theory. The Wagner effect which accounts for the change in the effective torsional stiffness due to the axial stresses is considered. 2. Based on the power series expansions of displacement components, the numerical method to exactly evaluate the element stiffness matrix of shear deformable thin-walled composite beam with arbitrary lay-up is presented. 3. For comparison, the finite element procedure based on the Largangian interpolation polynomial is presented. 4. To demonstrate accuracy and validity of this study, the critical buckling loads of I-beams with symmetric and arbitrary layups are evaluated and compared with the results from available references and the finite element solutions using isoparametric beam elements and shell elements from ABAQUS [24]. 2. General theory for shear deformable composite beam 2.1. Variational formulation In this study, two coordinate systems are used as shown in Fig. 1. The first coordinate system is the local coordinate (g, s, x), wherein the g axis is normal to the middle surface of the element and the s axis is tangent to the middle surface and is directed along the contour line of the cross-section. The second coordinate system is the orthogonal Cartesian coordinate (x, y, z), for which y and z axes lie in the plane of the cross-section and the x axis is parallel to the longitudinal axis of the beam; r is the g coordinate and q is the s coordinate of Q; x1 is the rotation of the cross-section about the axis of the beam and w is the angle tangent to the contour makes with y axis. To derive the composite beam theory, the following assumptions are made: 1. The strains are assumed to be small. 2. The beam is linearly elastic and prismatic. 3. The cross-section is assumed to maintain its shape during deformation, so that there is no distortion.
s, v
w Q q
,u
r
P O
Fig. 1. Pictorial definitions of coordinates in thin-walled cross-section.
102
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
4. The transverse shear strains and the warping shear are incorporated. It is assumed that they are uniform over the crosssections. From the study by Kim [25], the displacement fields of any generic point Q on the profile section can be expressed in terms of the seven displacement parameters as follows:
u ¼ U y sin w U z cos w x1 q @u @s w ¼ U x x3 y þ x2 z þ f / gðx3 sin w þ x2 cos w þ fqÞ
v ¼ Uy cos w þ Uz sin w þ x1 r g
dPE ¼
Z
1 2
L
Z A
o
dPG ¼
L
h
PfU 0y dU 0y
þ
U 0z dU 0z
þ ðx
þ
U 0y d
Y
! dt ¼ 0
ð8Þ
ext
t2
Pdt ¼
Z
t2 t1
Z
o
L
h
F 1 dU 0x þ F 2 dðU 0y x3 Þ þ F 3 dðU 0z þ x2 Þ
þ M2 dx02 M 3 dx03 þ M/ df 0 þ Tdðf þ x01 Þ þ Mt dðx01 f Þ þ PfU 0y dU 0y þ U 0z dU 0z þ ðx01 dU 0y þ U 0y dx01 Þzp ðx01 dU 0z þ U 0z dx01 Þyp i þM p x01 dx01 Te dDTe dxdt ¼ 0
ð4Þ
where yp and zp are coordinates of the pole P as shown in Fig. 1 and the following geometric relations are used.
ð9Þ
The constitutive equations [27] of the kth lamina in the local coordinate system of the cross-section are written as
rx rsx
"
k ¼
Q 11
Q 16
Q 16
Q 66
#k
ex csx
ð10Þ
where Q ij are the transformed reduced stiffness coefficients and include the material properties of each lamina. The stress–strain relationship can be simplified by adopting the zero hoop stress assumption (rs = 0) or the zero hoop strain one (es = 0) and the detailed explanation can be found in Ref. [27]. 2.2. Stability equations and force–displacement relationships The stability equations of the present study can be obtained by integrating the derivatives of the varied quantities by parts in Eq. (9), and collecting the coefficients of dUx, dUy, dUz, dx1, dx2, dx3, and df.
E11 U 00x þ E16 U 00y þ E17 U 00z þ ðE15 þ E18 Þx001 þ E12 x002 þ E17 x02 E13 x003 E16 x03 þ E14 f 00 þ ðE18 E15 Þf 0 ¼ 0
ð11aÞ
E16 U 00x þ E66 U 00y þ E67 U 00z þ ðE56 þ E68 Þx001 þ E26 x002 þ E67 x02 E36 x003 E66 x03 þ E46 f 00 þ ðE68 E56 Þf 0 þ PU 00y þ Pzp x001 ¼ 0
0 1 Þzp
x
o ðx01 dU 0z þ U 0z dx01 Þyp g þ Mp x01 dx01 dx
Z
ð3Þ
0 0 1 dU y
G
ð1cÞ
where A is the area of the cross-section. Variation of the potential energy can be expressed by substituting the assumed displacement fields in Eqs. (1a) and (1b) into Eq. (3) as follows:
Z
E
t1
ð2Þ
P 02 ðu þ v 02 ÞdAdx A
þ
Y
d
where L is the length of beam; F1 is the axial force; F2 and F3 are the shear forces in y and z directions, respectively; M2 and M3 are the bending moments about y and z axes, respectively; M/ is the bimoment; T and Mt are the two contributions to the total twisting moment; eox is the axial strain; jy and jz are the biaxial curvatures in y and z directions, respectively; j/ and jsx are the warping curvature and twisting curvature, respectively; coyx and cozx are the transverse shear strains; co/ is the warping shear strain. Now the potential energy due to the initial axial force P acting at the centroid is given as follows:
PG ¼
d t1
Y
ð1bÞ
ðF 1 deox þ F 2 dcoyx þ F 3 dcozx þ M2 djy þ M3 djz þ M/ dj/
þ Tdco/ þ M t djsx Þdx
t2
Substituting Eqs. (2), (4), and (7) into the extended Hamilton’s principle in Eq. (8), the following weak statement is obtained.
L
o
Z
ð1aÞ
where Ux, Uy, and Uz are the rigid body displacements; x1 is the rotation of cross-section; x2 and x3 are the rotations about y and z axes, respectively; f is the warping parameter. The variation of strain energy for the thin-walled composite beams considering the transverse shear and the warping induced shear deformation is expressed as follows [25]:
Z
In order to derive the stability equations and the force–displacement relationships of the laminated composite beam, the extended Hamilton’s principle is used by
ð11bÞ
E17 U 00x þ E67 U 00y þ E77 U 00z þ ðE57 þ E78 Þx001 þ E27 x002 þ E77 x02 E37 x003 E67 x03 þ E47 f 00 þ ðE78 E57 Þf 0 þ PU 00z Pyp x001 ¼ 0
ð11cÞ
ðE15 þ E18 ÞU 00x þ ðE56 þ E68 ÞU 00y þ ðE57 þ E78 ÞU 00z þ ð2E58 þ E55 þ E88 Þx001
y yp ¼ q cos w þ r sin w
ð5aÞ
þ ðE25 þ E28 Þx002 þ ðE57 þ E78 Þx02 ðE35 þ E38 Þx003 ðE56 þ E68 Þx03
z zp ¼ q sin w r cos w
ð5bÞ
þ ðE45 þ E48 Þf 00 ðE55 E88 Þf 0 þ Pzp U 00y Pyp U 00z þ Mp x001 ¼ 0
In Eq. (4), the stress resultant Mp which is known as the Wagner effect is evaluated by
Mp ¼
Z A
P 2 ðq þ r 2 þ g2 þ 2grÞdsdg A
ð6Þ
This property accounts for the change in the effective torsional stiffness due to the bending compressive and tensile stresses, which create a torque as the beam twists during buckling. For the isotropic beams, Eq. (4) is the same as the one in Kim et al. [26]. The variation of the potential energy due to external loads can be written as follows:
dPext ¼
dDTe Te
ð7Þ
where De and Te signify the member end displacement and force vector, respectively.
ð11dÞ
E12 U 00x E17 U 0x þ E26 U 00y E67 U 0y þ E27 U 00z E77 U 0z þ ðE25 þ E28 Þx001 ðE57 þ E78 Þx01 þ E22 x002 E77 x2 E23 x003 þ ðE37 E26 Þx03 þ E67 x3 þ E24 f 00 ðE25 E28 þ E47 Þf 0 þ ðE57 E78 Þf ¼ 0
ð11eÞ
E13 U 00x E16 U 0x þ E36 U 00y E66 U 0y þ E37 U 00z E67 U 0z þ ðE35 þ E38 Þx001 ðE56 þ E68 Þx01 þ E23 x002 þ ðE37 E26 Þx02 E67 x2 E33 x003 þ E66 x3 þ E34 f 00 ðE35 E38 þ E46 Þf 0 þ ðE56 E68 Þf ¼ 0 ð11fÞ E14 U 00x þ ðE15 E18 ÞU 0x þ E46 U 00y þ ðE56 E68 ÞU 0y þ E47 U 00z þ ðE57 E78 ÞU 0z þ ðE45 þ E48 Þx001 þ ðE55 E88 Þx01 þ E24 x002 þ ðE25 E28 þ E47 Þx02 þ ðE57 E78 Þx2 E34 x003 ðE35 E38 þ E46 Þx03 þ ðE68 E56 Þx3 þ E44 f 00 ðE55 2E58 þ E88 Þf ¼ 0
ð11gÞ
103
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
where Eij are the laminate stiffnesses which depend on the crosssection of the composite and given in Ref. [25]. The seven force– displacement relationships are as follows:
The solutions of seven displacement measures are taken as the following infinite power series.
F 1 ¼ E11 U 0x þ E16 U 0y þ E17 U 0z þ ðE15 þ E18 Þx01 þ E17 x2
Ux ¼
þ E12 x02 E16 x3 E13 x03 þ ðE18 E15 Þf þ E14 f 0
ð12aÞ
ð12bÞ
F 3 ¼ E17 U 0x þ E67 U 0y þ ðE77 þ PÞU 0y þ ðE78 þ E57 Pyp Þx01 ð12cÞ
M 1 ¼ ðE15 þ
Pyp ÞU 0z
þ ðE68 þ E56 þ
Pzp ÞU 0y
þ ðE78 þ E57
þ E78 Þx2 þ ðE25 þ E28 Þx02 ðE56 þ E68 Þx3 ðE35 þ E38 Þx þ ðE88 E55 Þf þ ðE45 þ E48 Þf M2 ¼
E12 U 0x
E26 U 0y
þ
þ
E27 U 0z
0
þ ðE25 þ E28 Þx þ E27 x2 ð12eÞ
E23 x02 þ E36 x3 þ E33 x03 þ ðE35 E38 Þf E34 f 0
ð12fÞ
þ E24 x02 E46 x3 E34 x03 þ ðE48 E45 Þf þ E44 f 0
ð12gÞ
The Eq. (11) is the most general form of the coupled buckling for the thin-walled laminated composite beam. If the cross-section is bisymmetric and the material is isotropic, in which the coupling effects are neglected, Eq. (11) can be simplified to the following differential equations [28].
EAU 00x ¼ 0
ð13aÞ
x
þ
¼0
ð13bÞ
GA3 ðU 00z þ x02 Þ þ PU 00z ¼ 0
ð13cÞ
GJ 001 þ GAr ð 001 þ f 0 Þ Mp EI2 002 GA3 ðU 0z þ 2 Þ ¼ 0 EI3 003 þ GA2 ðU 0y 3 Þ ¼ 0
x x x
x
1 X g n xn n¼0
Substituting Eq. (15) into Eq. (11) and shifting the index of power of xn, the stability equations can be compactly expressed in a matrix form as follows:
¼ Zn fan ; anþ1 ; bn ; bnþ1 ; cn ; cnþ1 ; dn ; dnþ1 ; en ; enþ1 ; fn ; fnþ1 ; g n ; g nþ1 gT ð16Þ
ð17Þ
From Eq. (16), the following relation is obtained.
00 1
x ¼0
ð13dÞ
x x
EI/ f 00 GAr ðx01 þ f Þ ¼ 0
¼ Zi Ni fai ; aiþ1 ; bi ; biþ1 ; ci ; ciþ1 ; di ; diþ1 ; ei ; eiþ1 ; fi ; fiþ1 ; g i ; g iþ1 gT ð18Þ
M / ¼ E14 U 0x þ E46 U 0y þ E47 U 0z þ ðE45 þ E48 Þx01 þ E47 x2
PU 00y
n¼0
f ¼
faiþ2 ; biþ2 ; ciþ2 ; diþ2 ; eiþ2 ; fiþ2 ; g iþ2 gT
0 1
M 3 ¼ E13 U 0x E36 U 0y E37 U 0z ðE35 þ E38 Þx01 E37 x2
0 3Þ
n¼0
1 X f n xn ;
a ¼ fao ; a1 ; bo ; b1 ; co ; c1 ; do ; d1 ; eo ; e1 ; fo ; f1 ; g o ; g 1 gT ð12dÞ
þ E22 x02 E26 x3 E23 x03 þ ðE28 E25 Þf þ E24 f 0
GA2 ðU 00y
x3 ¼
where Zn is the 7 14 matrix function. We then define the initial integration constant vector a as
þ ðE55 þ 2E58 þ E88 þ M p Þx01 þ ðE57
0 3
n¼0
1 X x2 ¼ en xn ;
fanþ2 ; bnþ2 ; cnþ2 ; dnþ2 ; enþ2 ; fnþ2 ; g nþ2 gT
þ E77 x2 þ E27 x02 E67 x3 E37 x03 þ ðE78 E57 Þf
E18 ÞU 0x
dn xn ;
n¼0
ð15Þ
þ E67 x2 þ E26 x02 E66 x3 E36 x03 þ ðE68 E56 Þf
þ E47 f 0
n¼0 1 X n¼0
F 2 ¼ E16 U 0x þ ðE66 þ PÞU 0y þ E67 U 0z þ ðE56 þ E68 þ Pzp Þx01 þ E46 f 0
x1 ¼
1 1 1 X X X an xn ; U y ¼ b n xn ; U z ¼ c n xn
The terms for ai+2, bi+2, ci+2, di+2, ei+2, fi+2, and gi+2 converge to zero as i ? 1. The displacement state vector in Eq. (14) is expressed with respect to a by using Eqs. (15) and (18) as follows:
d ¼ Xn a
ð19Þ
where Xn denotes the 14 14 matrix function with the coefficients of Ux, Uy, Uz, x1, x2, x3, and f. In each of these 14 solution sets, the calculation of the coefficients by the recursive relations in Eq. (18) is continued until the contribution of the next coefficient is less than an arbitrarily chosen small number. In this study, above symbolic calculations are performed with the help of the technical computer software Methematica [29]. Next, the initial integration constant vector a can be expressed with respect to 14 nodal displacement components at two ends of the laminated beam. For this, we consider the nodal displacement vector at a and b, which mean the two ends of member, is defined by
ð13eÞ
Ue ¼ hUa ; Ub iT
ð13fÞ
Substituting coordinates of the two ends of member (x = 0, L) into Eq. (19), the nodal displacement vector Ue is expressed as follows:
ð13gÞ
ð20Þ
where EA is the axial rigidity; GA2, and GA3 are the flexural shear rigidities with respect to y and z axes, respectively; GJ and GAr are the torsional rigidity and the warping shear rigidity, respectively; EI2 and EI3 are the flexural rigidities with respect to y and z axes, respectively; EI/ are the warping rigidity.
The displacement state vector consisting of 14 displacement components can be rewritten by using Eqs. (19) and (21) as follows:
3. Exact stiffness matrix of laminated composite beam
d ¼ Xn H1 Ue
3.1. Evaluation of exact displacement function In order to construct the member stiffness matrix of the laminated composite beam, the exact displacement functions are rigorously evaluated based on the seven stability equations derived in previous Chapter. For this, the displacement state vector consisting of 14 displacement parameters is considered as follows:
d¼
D
U x ; U 0x ; U y ; U 0y ; U z ; U 0z ;
0 1;
0 2;
0 0 3; f ; f
x1 ; x x2 ; x x3 ; x
ET
ð14Þ
Ue ¼ Ha
ð21Þ
ð22Þ
It should be noted that XnH1 in Eq. (22) is not an approximate displacement function but an exact one since the displacement state vector satisfies the homogenous form of the coupled equilibrium equations in Eq. (11) exactly. 3.2. Calculation of stiffness matrix In order to calculate the stiffness matrix of the composite member using the displacement function derived in this study, we consider the nodal force vector at two ends a and b defined by
104
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
Fe ¼ hFa ; Fb iT
ð23Þ
The force–deformation relations in Eq. (12) can be rewritten in a matrix form.
f ¼ Sd
ð24Þ
where S is the 7 14 matrix. Substitution of the displacement function in Eq. (22) into Eq. (24) leads to
f ¼ SXn H1 Ue
ð25Þ
Finally, the stiffness matrix of the fully coupled laminated composite beam is evaluated as
Fe ¼ KUe
ð26Þ
where
" K¼
SXn ð0ÞH1 SXn ðLÞH
# ð27Þ
1
The critical buckling loads Pcr of the beam are the values that cause the stiffness matrix for the element to become singular as Eq. (27). A search procedure is employed to find these values up to the desired accuracy. In this study, the Regular-Falsi method [30] is applied to ensure that none of the buckling loads is missed.
det jKj ¼ 0
ð28Þ
It is noteworthy that the stiffness matrix in Eq. (27) is formed by the shape functions which are exact solutions of the governing equations. Therefore, the thin-walled composite beam based on the stiffness matrix developed by this study eliminates discretization errors and is free from the shear and membrane locking. 4. Finite element formulation
The element displacement vector Uie and the force vector Fie for the isoparametric beam element are defined as
Uie ¼ ½U 1 ; U 2 n
U ¼
½U nx ; U ny ; U nz ;
ð30aÞ n 1;
n 2;
x x x
n n T 3; f ;
n ¼ 1; 2
ð30cÞ ð30dÞ
where Un and Fn are the nodal displacement and force vectors of the isoparametric beam element, respectively. Substituting the shape functions and the cross-sectional properties into Eq. (9) and integrating along the element length, the total potential energy of the finite composite beam element is obtained in a matrix form as
Y
¼
1 T U ðKie þ Kig ÞUie UTie Fie 2 ie
ð31Þ
where Kie and Kig are the element elastic and geometric stiffness matrices in local coordinate, respectively. When equal interpolation of all seven displacement parameters and full integration is used to evaluate the stiffness matrices Kie and Kig, the beam element is known to behavior very stiff. This behavior is due to the instability of the element to represent a constant state of the transverse shear strains coyx ; cozx and the warping shear strain co/ , and the phenomenon is known as the shear locking. To avoid locking, the shear strains are approximated, while using equal interpolation of the seven displacement parameters as a constant by using the reduced integration of the shear stiffnesses. Therefore, in the present work, the reduced integration scheme in which the quadrature rule is taken to be one order lower than the normal one is employed. For the entire structure, the assembly of element stiffness matrix leads to the equilibrium matrix equation in a global coordinate system based on the coordinate transformation.
For comparison, the finite element formulation based on the shear deformable isoparametric beam element with seven DOFs per a node is presented. The present finite beam element uses the same shape functions for all translational, rotational, and warping displacements. The displacements can be expressed as follows:
l ð1 þ r Þ 2 2 X Ui ¼ Rn ðr ÞU ni
x¼
ð30bÞ
Fie ¼ ½F 1 ; F 2 h iT Fn ¼ F n1 ; F n2 ; F n3 ; M n1 ; Mn2 ; M n3 ; Mn/
1
ð29aÞ i ¼ x; y; z
ð29bÞ
j ¼ 1; 2; 3
ð29cÞ
n¼1
xj ¼
2 X
Rn ðr Þxnj
n¼1
f ¼
2 X Rn ðr Þf n
ð29dÞ
n¼1
where U ni ; xni , and fn are the translational and rotational displacements and the warping parameter at node n, respectively; Rn is the Lagrangian interpolation function whose detailed expression is presented in Ref. [31]; r⁄ is a natural coordinate that varies from 1 to +1.
2
Fig. 2. Shape of cross-section under consideration.
Table 1 Critical buckling loads of the SS bisymmetric I-beams (106 N), (L = 6 m). Stacking sequence
[0]4 [45/ 45]S [0/90]S
Machado and Cortı´nez [32]
This study Without shear
With shear
rs = 0
rs = 0
42.69 4.45 22.90
33.18 4.44 19.85
Vo and Lee [11]
Without shear
With shear
es = 0
rs = 0
rs = 0
rs = 0
33.34 12.99 19.87
42.11 4.45 22.57
33.18 4.44 19.84
30.38 4.41 20.63
105
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
ðKE þ kKG ÞU ¼ F ¼ 0
ð32Þ
Table 2 Critical buckling loads of the CF bisymmetric I-beams (N), (L/h = 20).
where KE and KG are the global elastic and geometric stiffness matrices, respectively; F is the global nodal force vector; k is the load parameter under the assumption of proportional loading.
Stacking sequence
[0]16 [15/15]4S [30/30]4S [45/45]4S [60/60]4S [75/75]4S [0/90]4S
5. Numerical examples In order to illustrate accuracy and validity of the current study, numerical analyses have been performed for the coupled stability analysis of composite I-beams with the bi-symmetric and monosymmetric cross-sections. The terms ‘‘bisymmetric’’ and ‘‘monosymmetric’’ are defined in the geometry of the cross-section not the material point of view. The solutions obtained from this study are compared with those from other researchers and the finite element solutions using the isoparametric beam elements and the shell elements by ABAQUS.
Without shear rs = 0
With shear
rs = 0
es = 0
5755.2 5199.8 3861.0 2672.7 2114.7 1948.3 3857.8
5737.5 5187.7 3856.0 2670.6 2113.2 1946.6 3849.8
5859.1 5413.1 4325.2 3160.2 2370.6 2031.3 3910.8
Vo and Lee [11]
ABAQUS
5741.5 5189.0 3854.5 2668.4 2111.3 1945.1 3829.8
5719.2 5175.3 3853.3 2671.4 2113.4 1945.8 3843.4
5.1. Bisymmetric I-beams with symmetric lay-up The example considered is the simply supported (SS) bisymmetric I-beam, as shown in Fig. 2, whose geometric properties
0.5
7
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
0.4
6
Dimensionless buckling load, Pcr*
Dimensionless buckling load, Pcr*
This study
0.3
0.2
0.1
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
5
4
3
2
1
0.0
0 0
15
30
45
60
75
90
0
30
45
60
Fiber angle, (degree)
(a) CF beams
(c) CS beams
75
90
7
7
6
5
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
6
Without shear (L/h=5) With shear (L/h=5) Without shear (L/h=20) With shear (L/h=20)
Dimensionless buckling load, Pcr*
Dimensionless buckling load, Pcr*
15
Fiber angle, (degree)
4
3
2
5
4
3
2
1
1
0
0 0
15
30
45
60
75
90
0
15
30
45
60
Fiber angle, (degree)
Fiber angle, (degree)
(b) SS beams
(d) CC beams
Fig. 3. Critical buckling loads of the bisymmetric I-beams with respect to the fiber angle change.
75
90
106
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
are as follows: L = 6000 mm, b1 = b2 = 300 mm, h = 600 mm, and t = 30 mm. The cross-section is made with graphite-epoxy (AS4/ 3501) material and its material constants are as follows: E1 = 144 GPa, E2 = E3 = 9.65 GPa, G12 = G13 = 4.14 GPa, G23 = 3.45 GPa, m12 = m13 = 0.3, m23 = 0.5. In which subscripts ‘1’ and ‘2’ correspond to directions parallel and perpendicular to fibers, respectively. All constituent flanges and web are made of four plies and are assumed to be symmetrically laminated with respect to its middle plane. The critical buckling loads obtained from the present analysis using a single element are compared with the finite element results of Vo and Lee [11] and Machado and Cortı´nez [32] in Table 1. It is observed that the present solutions applying the zero hoop stress assumption (rs = 0) are in good agreement with results in Refs. [11,32]. Next, the bisymmetric I-beams made with glass–epoxy material with the following elastic constants: E1 = 53.78 GPa, E2 = E3 =
17.93 GPa, G12 = G13 = 8.96 GPa, G23 = 3.45 GPa, m12 = m13 = 0.25, m23 = 0.34 is considered. The total 16 layers with equal thickness are considered in two flanges and web and the sectional properties are as follows: b1 = b2 = 25 mm, h = 50 mm, and t = 2.08 mm. For clamped-free (CF) beams with L/h = 20, the buckling loads from this study with and without shear deformation effect are presented in Table 2. For comparison, the finite element solutions by Vo and Lee [11] and ABAQUS are given in Table 2. For ABAQUS calculation, a total of 240 nine node shell elements (S9R5) (six through the cross-section) are used to obtain the results. It can be found from Table 2 that the correlation of the present results considering the shear based on rs = 0 assumption and the solutions from ABAQUS and Vo and Lee [11] is seen to be excellent for all lay-ups considered. It can also be found that the analysis based on es = 0 assumption is seen to overestimate the buckling loads by a maximum 18.3% at w = 45°.
100
1.6
90
70
Dimensionless buckling load, Pcr*
Effect of shear, (%)
1.4
CF SS CS CC
80
60 50 40 30 20
E1/E2 1 3 5 10
1.2 1.0 0.8 0.6 0.4 0.2
10
0.0
0 0
15
30
45
60
75
90
0
15
30
45
60
Fiber angle, (degree)
Fiber angle, (degree)
(a) L h = 5
(a) Critical buckling load
75
90
75
90
0.20
6
0.19
Effect of shear, (%)
4
Dimensionless buckling load, Pcr*
CF SS CS CC
5
3
2
1
E1/E2 1 3 5 10
0.18 0.17 0.16 0.15 0.14 0.13 0.12
0 0
15
30
45
60
75
90
0
15
30
45
60
Fiber angle, (degree)
Fiber angle, (degree)
(b) L h = 20
(b) Region of interest
Fig. 4. Effect of shear for the bisymmetric I-beams with respect to the fiber angle change.
Fig. 5. Critical buckling loads for the flexural mode in y axis of the CF bisymmetric Ibeams with various modulus ratios, (L/h = 15).
107
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
To investigate the effects of boundary condition and shear deformation on the critical buckling loads of the laminated composite beam, the dimensionless buckling loads for four boundary conditions are plotted in Fig. 3 with respect to the fiber angle change. The dimensionless buckling load parameter 2 P cr ¼ P cr L2 =ðAE2 h Þ is introduced for simplicity and the boundary conditions under consideration are: the clamped–free (CF), simply–simply (SS), clamped–simply (CS), and clamped–clamped (CC) boundary conditions. It is seen that the buckling loads of beams without shear are the highest at w = 0°, both of short and long beams for all boundary conditions. Whereas for CC beam with shear, it is the highest around w = 14° for L/h = 5 beam with shear. Fig. 4 shows the variation of the effect of shear on the buckling loads for L/h = 5 and 20 beams. It is observed that the shear effects are the largest at w = 0° and the smallest at w = 50° for both of beams regardless of the boundary conditions. Moreover the maximum shear effect, as can be seen in Fig. 4a, is 5% for CC beam with L/h = 20 which may be relatively slender. Therefore the shear effect 5.0
Dimensionless buckling load, Pcr*
4.5 E1/E2
4.0
1 3 5 10
3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0
15
30
45
60
75
90
Fiber angle, (degree)
(a) Critical buckling load
should be included in the flexural–torsional stability analysis of the laminated composite beam. The effects of the modulus ratio E1/E2 on the critical buckling loads of the CF beams with L/h = 15 are presented in Figs. 5–7. It is seen that the maximum buckling loads of the isotropic beam (E1/E2 = 1) are obtained at w = 45° for the flexural modes in y and z axes as shown in Figs. 5 and 6. While those are exhibited at w = 0° for the anisotropic beams. It is also seen that the critical buckling load increases as the modulus ratio increases in fiber angles less than 60°. Whereas it is reverse in region beyond w = 60°. For the case of torsional mode as shown in Fig. 7, the torsional buckling loads are maximum and minimum at w = 90° for the isotropic and anisotropic beams, respectively. Fig. 8 shows the effect of shear on the flexural and torsional buckling loads. It is evident that the shear effect on the flexural buckling mode about the strong axis is larger than that about the weak axis. It is interesting to observe that for the isotropic beam, the maximum shear effect takes place at w = 45° for the two flexural modes and torsional mode. On the other hand, for anisotropic beams, it occurs at w = 0° for all buckling modes. It is also observed that the shear effect increases as E1/E2 increases in fiber angles w 6 25° and it decreases with increase of E1/E2 in w > 25°. Finally, the convergence study for the isoparametric beam element using the Lagrangian interpolation polynomials and the beam element developed by the present stiffness matrix method is performed and given in Table 3. In case of the stiffness matrix method, the system matrix can be obtained by the usual assembling process from the conditions of equilibrium and compatibility at the nodes. It can be found from Table 3 that the solutions developed by this study using one element are the same as results using three elements, and are in an excellent agreement with those using at least 60 isoparametric beam elements. It should be noted that both the present stiffness matrix method and the isoparametric beam formulation use the same total potential energy. The difference of two methods is that the current element is based on the shape functions which satisfy the homogeneous forms of the equilibrium equations exactly. Therefore it is possible to obtain the solutions though only a minimum number of beam elements are used. However, the solutions obtained from the isoparametric beam elements using the Lagrangian interpolation function which satisfies only displacement continuity at nodal point are approximate. Resultantly, as a number of isoparametric beam element
0.75 2.0 1.8
E1/E2 1 3 5 10
0.65
Dimensionless buckling load, Pcr*
Dimensionless buckling load, Pcr*
0.70
0.60
0.55
0.50
0.45
E1/E2 1 3 5 10
1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2
0.40 0
15
30
45
60
75
90
Fiber angle, (degree)
(b) Region of interest Fig. 6. Critical buckling loads for the flexural mode in z axis of the CF bisymmetric Ibeams with various modulus ratios, (L/h = 15).
0.0 0
15
30
45
60
75
90
Fiber angle, (degree) Fig. 7. Critical buckling loads for the torsional mode of the CF bisymmetric I-beams with various modulus ratios, (L/h = 15).
108
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
2.0
14
1.8 12 1 3 5 10
1.4 1.2
1 3 5 10
10
Effect of shear, (%)
Effect of shear, (%)
E1/E2
E1/E2
1.6
1.0 0.8 0.6
8
6
4
0.4 2 0.2 0
0.0 0
15
30
45
60
75
90
0
15
30
45
60
75
Fiber angle, (degree)
Fiber angle, (degree)
(a) Flexural mode in y axis
(b) Flexural mode in z axis
90
2.0 1.8 E1/E2
Effect of shear, (%)
1.6
1 3 5 10
1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0
15
30
45
60
75
90
Fiber angle, (degree)
(c) Torsional mode Fig. 8. Effect of shear of the CF bisymmetric I-beams with various modulus ratios, (L/h = 15).
Table 3 Convergence of the finite beam element and the present solutions for the buckling loads of the CF bisymmetric I-beam with [0/45/90/45]2S lay-up (N), (L/h = 20). Mode
1 2 3
No. of isoparametric beam elements
No. of elements by this study
ABAQUS
5
10
30
60
1
2
3
3379.9 11732.6 11818.3
3338.2 11589.1 11776.7
3326.0 11547.1 11764.6
3324.9 11543.2 11763.5
3324.5 11541.9 11763.1
3324.5 11541.9 11763.1
3324.5 11541.9 11763.1
used in the beam structure increases, its solutions converge to those from the present study using a minimum number of beam elements.
5.2. Monosymmetric I-beams with symmetric and arbitrary lay-ups In this example, the monosymmetric I-beams with symmetric and arbitrary lay-ups are considered. The material is the glass–
3321.3 11350.0 11478.0
epoxy which is used in Section 5.1 and the geometric properties are b1 = 15 mm, b2 = 25 mm, h = 50 mm, and t = 2.08 mm. First, for the CF monosymmetric I-beam with symmetric lay-up, the coupled buckling loads are compared with those from ABAQUS analysis in Table 4 for various stacking sequences. In Table 4, the term ‘‘NW’’ means the buckling load from neglecting warping. It can be found from Table 4 that the good correlation of the current study including rs = 0 assumption and considering warping effect with the ABAQUS analysis is obtained over the whole range of fiber
109
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111 Table 4 Critical buckling loads of the CF monosymmetric I-beams (N), (L/h = 20). Stacking sequence
This study
[0]16 [15/15]4S [30/30]4S [45/45]4S [60/60]4S [75/75]4S [0/90]4S
1.2 ABAQUS
Without shear rs = 0
With shear
rs = 0
rs = 0(NW)
es = 0
3005.9 2854.6 2201.0 1562.4 1241.7 1134.6 2114.0
2992.5 2809.4 2199.0 1561.5 1240.9 1133.8 2111.1
3492.7 3157.6 2346.5 1625.0 1285.8 1184.4 2343.2
3046.4 2917.2 2447.5 1834.9 1386.5 1181.2 2141.0
1.0 0.8
2812.7 2664.9 2129.2 1534.8 1226.6 1122.4 2055.7
0.6
U U U
ω ω ω f
0.4 0.2 0.0
Table 5 Critical buckling loads of the CF monosymmetric I-beam with [0/30/60/90] lay-up (N), (L/h = 15). Mode
1 2 3
This study
ABAQUS
Without shear
With shear
rs = 0
rs = 0
es = 0
3171.1 10951.8 17321.8
3135.9 10900.1 16851.8
3271.9 11176.9 17688.9
-0.2 -0.4 0.00
0.75
1.00
0.75
1.00
(a) Mode 1 2921.1 11189.0 16391.0
Table 6 Critical buckling loads of the CC monosymmetric I-beam with [0/30/60/90] lay-up (N), (L/h = 15).
1 2 3
0.50
x/ L
1.2 1.0
Mode
0.25
0.8 0.6
Without shear
With shear
rs = 0
rs = 0
es = 0
29792.7 57543.3 93819.4
28842.0 54631.8 89936.0
30017.1 56952.8 94770.5
U U U
ω ω ω f
0.4 0.2 0.0
angle. The maximum deviation with es = 0 assumption is about 17.5% at 45° of fiber angle and the effect of warping is the largest at w = 0° with 16.7% and the smallest at w = 60° with 3.6%. It can also be noticed that the present approach exhibits better results than that of classical beam assumption which does not consider shear effect. The purpose of our final example is to evaluate exactly the coupled buckling load of the monosymmetric I-beam with arbitrary lay-up. We consider the monosymmetric I-beam with the same dimensions used in the previous example, but four layers of [0/ 30/60/90] with equal ply thicknesses are used for both flanges and web. Since the lay-up is nonsymmetric with respect to the middle plane of the plate wall, the extension-shear (A16 and A26), bending–twisting (D16 and D26), and extension–twisting and bending–shearing (B16 and B26) couplings have non-zero values. Resultantly all coupling laminate stiffnesses in Eqs. (11) and (12) are non-zero except for E12,E67, and E78. For CF and CC beams with L/ h = 15, the lowest three coupled buckling loads are presented in Tables 5 and 6, respectively. It is seen that the results from the present rs = 0 assumption are in close with ABAQUS solutions for CF beam. The shear effect reduces the coupled buckling loa ds about 2.8% in the 3rd mode for CF beam and 5.3% in the 2nd mode for CC beam. It should be noted that the present solutions obtained from one element are exact for the higher buckling modes as well as the lower ones, while the considerable numbers of finite beam elements based on approximate shape functions may be required to achieve sufficiently accurate results in the higher modes. Furthermore for beams with L/h = 10, the buckling mode shapes corresponding to the lowest two buckling loads are illustrated in
-0.2 -0.4 0.00
0.25
0.50
x/ L
(b) Mode 2 Fig. 9. Coupled buckling mode shapes of CF monosymmetric I-beam with [0/30/60/ 90] lay-up, (L/h = 10).
Figs. 9 and 10 for CF and CC beams, respectively. From Figs. 9 and 10, the flexural–torsional coupled buckling modes and the flexural–torsional-shearing coupled ones can be observed for CF and CC beams, respectively.
6. Conclusion The super convergent finite composite beam elements considering shear effects has been developed for the coupled stability analysis. The stability equations and force–displacement relationships are derived by using the extended Hamilton’s principle and the element stiffness matrix is rigorously evaluated based on the power series expansions of displacement components. For comparison, the finite element model using the Lagrangian interpolation polynomial is presented. An advantage, which is often overlooked but may be more important, is that the present method can provide
110
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111
1.2 U U U
1.0
ω ω ω
0.8 0.6
f
0.4 0.2 0.0 -0.2 -0.4 -0.6 0.00
0.25
0.50
0.75
1.00
x/ L
(a) Mode 1
(2) For both of the shorter and longer span beams with bisymmetric cross-section, the minimum shear effect is reached at w = 55° regardless of the boundary conditions. (3) The torsional buckling loads of beams with high modulus ratio are larger than those with small modulus ration whole fiber angles except w = 90°. But it is not always a case of the flexural buckling behavior. (4) For flexural and torsional bucking behavior, the shear effect with high modulus ratio is not always larger than that with small modulus ratio. (5) The present numerical procedure provides a refined method for not only the evaluation of the stiffness matrix of shear deformable laminated composite beam but also general solutions of simultaneous ordinary differential equations of the higher order. This thin-walled composite beam element eliminates discretization errors and is free from locking phenomenon since the displacement state vector satisfies the homogenous form of the simultaneous ordinary differential equations. Furthermore this study is capable of predicting an infinite number of buckling loads of composite beams by means of a finite number of coordinates.
1.2 1.0 0.8 0.6 0.4
U U U
Acknowledgements This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2012R1A1A2007054).
ω ω ω f
0.2
References 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 -1.2 0.00
0.25
0.50
0.75
1.00
x/ L
(b) Mode 2 Fig. 10. Coupled buckling mode shapes of CC monosymmetric I-beam with [0/30/ 60/90] lay-up, (L/h = 10).
benchmark results when compared with other results obtained by the finite element or other approximate methods. Through the numerical examples, the present study is validated comparing the coupled buckling responses with results from other researchers and the finite element solutions using the isoparametric beam elements and the shell elements of ABAQUS. Good correlation is achieved for various beams with different cross-sections, lamination schemes, boundary conditions considered in this paper. The conclusions drawn from this study are as follows: (1) The flexural buckling behavior is affected by the shear and the span-to-height ratio as the fiber angle varies: for shorter span beam with CC boundary condition and with shear effect, the maximum buckling load is obtained at other fiber angle except w = 0°. Whereas for longer span beams with and without shear, and for shorter beam without shear, the maximum buckling loads occur at w = 0°.
[1] Vlasov VZ. Thin walled elastic beams. 2nd ed. Jerusalem: Israel Program for Scientific Transactions; 1961. [2] Gjelsvik A. The theory of thin-walled bars. Wiley; 1981. [3] Bauld NR, Tzeng L. A Vlasov theory for fiber-reinforced beams with thin-walled open cross sections. Int J Solids Struct 1984;20:277–97. [4] Kollá LP. Flexural–torsional buckling of open section composite columns with shear deformation. Int J Solids Struct 2001;38:7525–41. [5] Davalos JF, Qiao P, Salim HA. Flexural–torsional buckling of pultruded fiber reinforced plastic composite I-beams: experimental and analytical evaluations. Compos Strcut 1997;38:241–50. [6] Qiao P, Zou G, Davalos JF. Flexural–torsional buckling of fiber-reinforced plastic composite cantilever I-beams. Compos Strcut 2003;60:205–17. [7] Shan L, Qiao P. Flexural–torsional buckling of fiber-reinforced plastic composite open channel beams. Compos Strcut 2005;68:211–24. [8] Sherbourne AN, Kabir MZ. Shear strain effects in lateral stability of thin-walled fibrous composite beams. J Eng Mech 1995;121:640–7. [9] Pandey MD, Kabir MZ, Sherbourne AN. Flexural–torsional stability of thinwalled composite I-section beams. Compos Eng 1995;5:321–42. [10] Cortı´nez V, Piovan MT. Vibration and buckling of composite thin-walled beams with shear deformability. J Sound Vib 2002;258:701–23. [11] Vo TP, Lee J. On sixfold coupled buckling of thin-walled composite beams. Compos Struct 2009;90:295–303. [12] Vo TP, Lee J. Flexural–torsional coupled vibration and buckling of thin-walled open section composite beams using shear-deformable beam theory. Int J Mech Sci 2009;51:631–41. [13] Back SY, Will KM. Shear-flexible thin-walled element for composite I-beams. Eng Struct 2008;30:1447–58. [14] Suresh R, Malhotra SK. Some studies on buckling of laminated composite thin walled box beams. Compos Struct 1998;40:267–75. [15] Lin ZM, Polyzois D, Shah A. Stability of thin-walled pultruded structural members by finite element method. Thin-Walled Struct 1996;24:1–18. [16] Cortı´nez V, Piovan MT. Stability of composite thin-walled beams with shear deformability. Comput Struct 2006;84:978–90. [17] Silva NMF, Silvestre N, Camotim D. GBT formulation to analyze the buckling behavior of FRP composite open-section thin-walled columns. Compos Strcut 2010;93:79–92. [18] Silva NMF, Silvestre N. On the influence of material couplings on the linear and buckling behavior of I-section composite columns. Int J Strcut Stab Dyn 2007;7:243–72. [19] Shield CK, Morey TA. Kinematic theory for buckling of open and closed section thin-walled composite beams. J Eng Mech 1997;123:1070–81. [20] Bhaskar K, Librescu L. Buckling under axial compression of thin-walled composite beams exhibiting extension-twist coupling. Compos Struct 1995;31:203–12.
N.-I. Kim, D.-H. Choi / Composites: Part B 44 (2013) 100–111 [21] Washizu K. Variational method in elasticity and plasticity. Pergammon Press; 1975. [22] Jun L, Hongxing H, Rongying S. Dynamic stiffness analysis for free vibrations of axially loaded laminated composite beams. Compos Struct 2008;84: 87–98. [23] Abramovich H, Eisenberger M, Shulepov O. Vibration and buckling of cross-ply nonsymmetric laminated composite beams. AIAA J 1996;34:1064–9. [24] ABAQUS. Standard user’s manual, Ver. 6.1. Hibbit, Kalsson & Sorensen Inc; 2003. [25] Kim NI. Shear deformable doubly- and mono-symmetric composite I-beams. Int J Mech Sci 2011;53:31–41.
111
[26] Kim MY, Chang SP, Kim SB. Spatial stability and free vibration of shear flexible thin-walled elastic beams. I: analytical approach. Int J Numer Meth Eng 1994;37:4097–115. [27] Jones RM. Mechanics of composite material. Taylor and Francis; 1999. [28] Kim MY, Kim NI, Yun HT. Exact dynamic and static stiffness matrices of shear deformable thin-walled beam-columns. J Sound Vib 2003;267:29–55. [29] Wolfram Mathematica 7. Wolfram Research Inc., IL; 2009. [30] Wendroff B. Theoretical numerical analysis. NY: Academic Press; 1966. [31] Bathe KJ. Finite element procedures. Englewood Cliffs, NJ: Prentice-Hall; 1996. [32] Machado SP, Cortı´nez VH. Non-linear model for stability of thin-walled composite beams with shear deformation. Thin-Walled Struct 2005;43:1615–45.